diff --git a/src/common/aka_static_if.hh b/src/common/aka_static_if.hh index cc77f848d..29a9bae34 100644 --- a/src/common/aka_static_if.hh +++ b/src/common/aka_static_if.hh @@ -1,29 +1,94 @@ -/** - * @file aka_static_if.hh - * - * @author Nicolas Richart - * - * @date creation: Sun Aug 13 2017 - * @date last modification: Wed Oct 25 2017 - * - * @brief if constexpr equivalent in c++14 - * - * @section LICENSE - * - * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) - * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) - * - * Akantu 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. - * - * Akantu 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 Akantu. If not, see . - * - */ +// Copyright (c) 2016 Vittorio Romeo +// License: AFL 3.0 | https://opensource.org/licenses/AFL-3.0 +// http://vittorioromeo.info | vittorio.romeo@outlook.com + +#ifndef __AKANTU_AKA_STATIC_IF_HH__ +#define __AKANTU_AKA_STATIC_IF_HH__ + +#include + +#define FWD(...) ::std::forward(__VA_ARGS__) + +namespace akantu { + +template auto static_if(TPredicate) noexcept; + +namespace impl { + template struct static_if_impl; + + template struct static_if_result; + + template auto make_static_if_result(TF && f) noexcept; + + template <> struct static_if_impl { + template auto & else_(TF &&) noexcept { + // Ignore `else_`, as the predicate is true. + return *this; + } + + template auto & else_if(TPredicate) noexcept { + // Ignore `else_if`, as the predicate is true. + return *this; + } + + template auto then(TF && f) noexcept { + // We found a matching branch, just make a result and + // ignore everything else. + return make_static_if_result(FWD(f)); + } + }; + + template <> struct static_if_impl { + template auto & then(TF &&) noexcept { + // Ignore `then`, as the predicate is false. + return *this; + } + + template auto else_(TF && f) noexcept { + // (Assuming that `else_` is after all `else_if` calls.) + + // We found a matching branch, just make a result and + // ignore everything else. + + return make_static_if_result(FWD(f)); + } + + template auto else_if(TPredicate) noexcept { + return static_if(TPredicate{}); + } + + template auto operator()(Ts &&...) noexcept { + // If there are no `else` branches, we must ignore calls + // to a failed `static_if` matching. + } + }; + + template + struct static_if_result : TFunctionToCall { + // Perfect-forward the function in the result instance. + template + explicit static_if_result(TFFwd && f) noexcept : TFunctionToCall(FWD(f)) {} + + // Ignore everything, we found a result. + template auto & then(TF &&) noexcept { return *this; } + + template auto & else_if(TPredicate) noexcept { + return *this; + } + + template auto & else_(TF &&) noexcept { return *this; } + }; + + template auto make_static_if_result(TF && f) noexcept { + return static_if_result{FWD(f)}; + } +} // namespace impl + +template auto static_if(TPredicate) noexcept { + return impl::static_if_impl{}; +} + +#undef FWD +} // namespace akantu + +#endif /* __AKANTU_AKA_STATIC_IF_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc index 139d64a96..a2030ce7f 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_exponential.cc @@ -1,349 +1,338 @@ /** * @file material_cohesive_exponential.cc * * @author Mauro Corrado * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Mon Jul 09 2012 * @date last modification: Tue Feb 20 2018 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_exponential.hh" #include "dof_synchronizer.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveExponential::MaterialCohesiveExponential( SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model, id) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable, "Beta parameter"); this->registerParam("exponential_penalty", exp_penalty, true, _pat_parsable, "Is contact penalty following the exponential law?"); this->registerParam( "contact_tangent", contact_tangent, Real(1.0), _pat_parsable, "Ratio of contact tangent over the initial exponential tangent"); // this->initInternalArray(delta_max, 1, _ek_cohesive); use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); if ((exp_penalty) && (contact_tangent != 1)) { contact_tangent = 1; AKANTU_DEBUG_WARNING("The parsed paramter is forced to " "1.0 when the contact penalty follows the exponential " "law"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); - auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); - auto normal_it = normal.begin(spatial_dimension); - auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); - auto delta_max_it = delta_max(el_type, ghost_type).begin(); - auto delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); /// compute scalars Real beta2 = beta * beta; /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++delta_max_it, ++delta_max_prev_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \beta^2 \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); if ((normal_opening_norm < 0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { Vector op_n(*normal_it); op_n *= normal_opening_norm; Vector delta_s(*opening_it); delta_s -= op_n; delta = tangential_opening_norm * beta; computeCoupledTraction(*traction_it, *normal_it, delta, delta_s, *delta_max_it, *delta_max_prev_it); computeCompressiveTraction(*traction_it, *normal_it, normal_opening_norm, *opening_it); } else computeCoupledTraction(*traction_it, *normal_it, delta, *opening_it, *delta_max_it, *delta_max_prev_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTraction( Vector & tract, const Vector & normal, Real delta, const Vector & opening, Real & delta_max_new, Real delta_max) { AKANTU_DEBUG_IN(); /// full damage case if (std::abs(delta) < Math::getTolerance()) { /// set traction to zero tract.clear(); } else { /// element not fully damaged /** * Compute traction loading @f$ \mathbf{T} = * e \sigma_c \frac{\delta}{\delta_c} e^{-\delta/ \delta_c}@f$ */ /** * Compute traction unloading @f$ \mathbf{T} = * \frac{t_{max}}{\delta_{max}} \delta @f$ */ Real beta2 = beta * beta; Real normal_open_norm = opening.dot(normal); Vector op_n_n(spatial_dimension); op_n_n = normal; op_n_n *= (1 - beta2); op_n_n *= normal_open_norm; tract = beta2 * opening; tract += op_n_n; delta_max_new = std::max(delta_max, delta); tract *= exp(1) * sigma_c * exp(-delta_max_new / delta_c) / delta_c; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressiveTraction( Vector & tract, const Vector & normal, Real delta_n, __attribute__((unused)) const Vector & opening) { Vector temp_tract(normal); if (exp_penalty) { temp_tract *= delta_n * exp(1) * sigma_c * exp(-delta_n / delta_c) / delta_c; } else { Real initial_tg = contact_tangent * exp(1) * sigma_c * delta_n / delta_c; temp_tract *= initial_tg; } tract += temp_tract; } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTangentTraction( const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array::matrix_iterator tangent_it = - tangent_matrix.begin(spatial_dimension, spatial_dimension); - - Array::matrix_iterator tangent_end = - tangent_matrix.end(spatial_dimension, spatial_dimension); - + auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); + auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = normal.begin(spatial_dimension); - auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); - auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); Real beta2 = beta * beta; /** * compute tangent matrix @f$ \frac{\partial \mathbf{t}} * {\partial \delta} = \hat{\mathbf{t}} \otimes * \frac{\partial (t/\delta)}{\partial \delta} * \frac{\hat{\mathbf{t}}}{\delta}+ \frac{t}{\delta} [ \beta^2 \mathbf{I} + * (1-\beta^2) (\mathbf{n} \otimes \mathbf{n})] @f$ **/ /** * In which @f$ * \frac{\partial(t/ \delta)}{\partial \delta} = * \left\{\begin{array} {l l} * -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad if * \delta \geq \delta_{max} \\ * 0 & \quad if \delta < \delta_{max}, \delta_n > 0 * \end{array}\right. @f$ **/ for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it) { Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); if ((normal_opening_norm < 0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { Vector op_n(*normal_it); op_n *= normal_opening_norm; Vector delta_s(*opening_it); delta_s -= op_n; delta = tangential_opening_norm * beta; computeCoupledTangent(*tangent_it, *normal_it, delta, delta_s, *delta_max_it); computeCompressivePenalty(*tangent_it, *normal_it, normal_opening_norm); } else computeCoupledTangent(*tangent_it, *normal_it, delta, *opening_it, *delta_max_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTangent( Matrix & tangent, const Vector & normal, Real delta, const Vector & opening, Real) { AKANTU_DEBUG_IN(); Real beta2 = beta * beta; Matrix J(spatial_dimension, spatial_dimension); J.eye(beta2); if (std::abs(delta) < Math::getTolerance()) { delta = Math::getTolerance(); } Real opening_normal; opening_normal = opening.dot(normal); Vector delta_e(normal); delta_e *= opening_normal; delta_e *= (1 - beta2); delta_e += (beta2 * opening); Real exponent = exp(1 - delta / delta_c) * sigma_c / delta_c; Matrix first_term(spatial_dimension, spatial_dimension); first_term.outerProduct(normal, normal); first_term *= (1 - beta2); first_term += J; Matrix second_term(spatial_dimension, spatial_dimension); second_term.outerProduct(delta_e, delta_e); second_term /= delta; second_term /= delta_c; Matrix diff(first_term); diff -= second_term; tangent = diff; tangent *= exponent; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressivePenalty( Matrix & tangent, const Vector & normal, Real delta_n) { if (!exp_penalty) delta_n = 0; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); Real normal_tg = contact_tangent * exp(1) * sigma_c * exp(-delta_n / delta_c) * (1 - delta_n / delta_c) / delta_c; n_outer_n *= normal_tg; tangent += n_outer_n; } INSTANTIATE_MATERIAL(cohesive_exponential, MaterialCohesiveExponential); } // akantu diff --git a/src/solver/sparse_matrix_inline_impl.cc b/src/solver/sparse_matrix_inline_impl.cc index 8a90f0860..81da7d999 100644 --- a/src/solver/sparse_matrix_inline_impl.cc +++ b/src/solver/sparse_matrix_inline_impl.cc @@ -1,36 +1,35 @@ /** * @file sparse_matrix_inline_impl.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Mon Jun 19 2017 * * @brief implementation of inline methods of the SparseMatrix class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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. * * Akantu 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 Akantu. If not, see . * */ - /* -------------------------------------------------------------------------- */ namespace akantu { inline void SparseMatrix::clearProfile() { this->nb_non_zero = 0; } } // akantu diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fbbf6db2c..d3c826166 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,67 +1,66 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Alejandro M. Aragón # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 -# @date last modification: Fri Jan 22 2016 +# @date last modification: Mon Feb 12 2018 # # @brief configuration for tests # # @section LICENSE # -# Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de -# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des -# Solides) +# Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # -# Akantu 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 +# Akantu 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. # -# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# Akantu 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 +# 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 +# You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_DIR} ) set(AKANTU_TESTS_FILES CACHE INTERNAL "") #=============================================================================== # List of tests #=============================================================================== add_akantu_test(test_common "Test the common part of Akantu") add_akantu_test(test_static_memory "Test static memory") add_akantu_test(test_fe_engine "Test finite element functionalties") add_akantu_test(test_mesh_utils "Test mesh utils") add_akantu_test(test_model "Test model objects") add_akantu_test(test_solver "Test solver function") add_akantu_test(test_io "Test the IO modules") add_akantu_test(test_contact "Test the contact part of Akantu") add_akantu_test(test_geometry "Test the geometry module of Akantu") add_akantu_test(test_synchronizer "Test synchronizers") add_akantu_test(test_python_interface "Test python interface") package_add_files_to_package( cmake/akantu_test_driver.sh cmake/AkantuTestsMacros.cmake ) package_is_activated(parallel _is_parallel) if (_is_parallel) option(AKANTU_TESTS_ALWAYS_USE_MPI "Defines if sequential tests should also use MPIEXEC" FALSE) mark_as_advanced(AKANTU_TESTS_ALWAYS_USE_MPI) endif() diff --git a/test/python_fe/dynamic_solver.py b/test/python_fe/dynamic_solver.py index 1d45d9eb5..42d9e532c 100644 --- a/test/python_fe/dynamic_solver.py +++ b/test/python_fe/dynamic_solver.py @@ -1,90 +1,101 @@ +# ------------------------------------------------------------------------------ +__author__ = "Nicolas Richart" +__copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ + " en Mécanique des Solides)" +__credits__ = ["Nicolas Richart"] +__license__ = "L-GPLv3" +__maintainer__ = "Nicolas Richart" +__email__ = "nicolas.richart@epfl.ch" +# ------------------------------------------------------------------------------ + import copy import numpy.linalg as npl import scipy.sparse as sp import scipy.sparse.linalg as spl from . import export from . import fe @export class DynamicSolver(fe.Solver): def __init__(self, model, **kwargs): opt = copy.copy(kwargs) self._delta_t = opt.pop("delta_t", 0.001) self._alpha = opt.pop("alpha", 1./2.) self._beta = opt.pop("beta", 1./2.) self._type = opt.pop("type", 'disp') self._tolerance = opt.pop("tolerance", 1e-10) self._max_nloops = opt.pop("max_iterations", 100) self._model = model self._J = sp.csr_matrix(self._model.K.shape) self._coeff = {'disp': {'disp': 1., 'velo': 1. / (self._alpha * self._delta_t), 'acce': 1. / (self._alpha * self._beta * self._delta_t ** 2)}, # NOQA: E501 'velo': {'disp': self._alpha * self._delta_t, 'velo': 1., 'acce': 1. / (self._beta * self._delta_t)}, 'acce': {'disp': self._alpha * self._beta * self._delta_t ** 2, # NOQA: E501 'velo': self._beta * self._delta_t, 'acce': 1.}} def _assembleResidual(self): self._r = self._model.f_ext - self._model.f_int - \ self._model.M * self._model.a C = self._model.C if C is not None: self._r -= C * self._model.v def _predictor(self): self._model.u += self._delta_t * self._model.v + \ self._delta_t ** 2 / 2. * self._model.a self._model.v += self._delta_t * self._model.a def _corrector(self, delta_): self._model.u += self._coeff[self._type]['disp'] * delta_ self._model.v += self._coeff[self._type]['velo'] * delta_ self._model.a += self._coeff[self._type]['acce'] * delta_ def _assembleJacobian(self): K = self._model.K e = self._coeff[self._type]['disp'] self._J = e * K C = self._model.C if C is not None: d = self._coeff[self._type]['velo'] self._J += d * C M = self._model.M if M is not None: c = self._coeff[self._type]['acce'] self._J += c * M self._model.applyDirichletBC() self._zero_rows(self._J, self._model.blocked) def solveStep(self): self._predictor() self._nloops = 0 converged = False while not converged: self._assembleJacobian() self._assembleResidual() delta_ = spl.spsolve(self._J, self._r) self._corrector(delta_) self._nloops += 1 error = npl.norm(delta_) converged = error < self._tolerance or \ self._nloops > self._max_nloops print("{0} {1} -> {2}".format(error, self._nloops, converged)) if self._nloops >= self._max_nloops: raise ValueError('The solver did not converge') @property def nloops(self): return self._nloops diff --git a/test/python_fe/fe.py b/test/python_fe/fe.py index 7420b4696..e80e3f872 100644 --- a/test/python_fe/fe.py +++ b/test/python_fe/fe.py @@ -1,49 +1,60 @@ +# ------------------------------------------------------------------------------ +__author__ = "Nicolas Richart" +__copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ + " en Mécanique des Solides)" +__credits__ = ["Nicolas Richart"] +__license__ = "L-GPLv3" +__maintainer__ = "Nicolas Richart" +__email__ = "nicolas.richart@epfl.ch" +# ------------------------------------------------------------------------------ + import abc class Solver(metaclass=abc.ABCMeta): @abc.abstractmethod def solveStep(self): pass def _zero_rows(self, A, rows): for r in rows: diag = A[r, r] # zeros the lines defined in b A.data[A.indptr[r]:A.indptr[r + 1]] = 0. A[r, r] = diag class Model(metaclass=abc.ABCMeta): @abc.abstractmethod def applyDirichletBC(self): pass @abc.abstractmethod def f_int(self): pass @property @abc.abstractmethod def f_ext(self): pass @property @abc.abstractmethod def K(self): pass @property @abc.abstractmethod def M(self): pass @property @abc.abstractmethod def C(self): pass @property @abc.abstractmethod def blocked(self): pass diff --git a/test/python_fe/static_solver.py b/test/python_fe/static_solver.py index 67c48be66..201e725e1 100644 --- a/test/python_fe/static_solver.py +++ b/test/python_fe/static_solver.py @@ -1,25 +1,36 @@ +# ------------------------------------------------------------------------------ +__author__ = "Nicolas Richart" +__copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ + " en Mécanique des Solides)" +__credits__ = ["Nicolas Richart"] +__license__ = "L-GPLv3" +__maintainer__ = "Nicolas Richart" +__email__ = "nicolas.richart@epfl.ch" +# ------------------------------------------------------------------------------ + import copy import scipy.sparse.linalg as spl from . import export from . import fe @export class StaticSolver(fe.Solver): def __init__(self, model): self._model = model def _assembleResidual(self): self._r = self._model.f_ext - self._model.f_int def _assembleJacobian(self): self._J = copy.copy(self._model.K) self._model.applyDirichletBC() self._zero_rows(self._J, self._model.blocked) def solveStep(self): self._assembleJacobian() self._assembleResidual() x = spl.spsolve(self._J, self._r) self._model.u += x diff --git a/test/python_fe/truss_fe.py b/test/python_fe/truss_fe.py index c37a90cf7..757a49c7a 100644 --- a/test/python_fe/truss_fe.py +++ b/test/python_fe/truss_fe.py @@ -1,121 +1,132 @@ +# ------------------------------------------------------------------------------ +__author__ = "Nicolas Richart" +__copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ + " en Mécanique des Solides)" +__credits__ = ["Nicolas Richart"] +__license__ = "L-GPLv3" +__maintainer__ = "Nicolas Richart" +__email__ = "nicolas.richart@epfl.ch" +# ------------------------------------------------------------------------------ + import numpy as np import copy import scipy.sparse as sp from . import export from . import fe @export class TrussFE(fe.Model): def __init__(self, **kwargs): opt = copy.copy(kwargs) self._A = opt.pop("A", 1.) self._E = opt.pop("E", 1.) self._L = opt.pop("L", 1.) self._rho = opt.pop("rho", 1.) self._F = opt.pop("F", {-10: [-1]}) self._blocked = opt.pop("blocked", ([0], [0])) self._Ne = opt.pop("Ne", 1) self._Le = self._L / self._Ne self._F_ext = np.zeros((self._Ne + 1)) self._F_int = np.zeros((self._Ne + 1)) self._u = np.zeros((self._Ne + 1)) self._v = np.zeros((self._Ne + 1)) self._a = np.zeros((self._Ne + 1)) self.assembleMass() def assembleMass(self): Me = np.array([[2., 1.], [1., 2.]]) Me *= self._rho * self._A * self._Le / 6 self._M = sp.lil_matrix((self._Ne + 1, self._Ne + 1)) for i in range(self._Ne): self._M[i:i+2, i:i+2] += Me self._M = self._M.tocsr() def assembleStiffness(self): Ke = np.array([[1., -1.], [-1., 1.]]) Ke *= self._E * self._A / self._Le self._K = sp.lil_matrix((self._Ne + 1, self._Ne + 1)) for i in range(self._Ne): self._K[i:i+2, i:i+2] += Ke self._K = self._K.tocsr() def applyNeumannBC(self): for force, nodes in self._F.items(): self._F_ext[nodes] = force def applyDirichletBC(self): for i, b in enumerate(self._blocked[0]): self._u[b] = self._blocked[1][i] def computeInternalForces(self): self._F_int = self._K * self._u @property def K(self): self.assembleStiffness() return self._K @property def M(self): return self._M @property def C(self): return None @property def u(self): self.applyDirichletBC() return self._u @u.setter def u(self, x): self._u = x self.applyDirichletBC() @property def v(self): return self._v @v.setter def v(self, v_): self._v = v_ @property def a(self): return self._a @a.setter def a(self, a_): self._a = a_ @property def f_int(self): self.computeInternalForces() return self._F_int @property def f_ext(self): self.applyNeumannBC() return self._F_ext @f_ext.setter def f_ext(self, f): self._F_ext = f @property def blocked(self): return self._blocked[0] @property def Le(self): return self._Le diff --git a/test/test_common/CMakeLists.txt b/test/test_common/CMakeLists.txt index 1a8bc7149..f95146d65 100644 --- a/test/test_common/CMakeLists.txt +++ b/test/test_common/CMakeLists.txt @@ -1,46 +1,37 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 -# @date last modification: Mon Dec 07 2015 +# @date last modification: Tue Dec 05 2017 # # @brief configurations for common tests # # @section LICENSE # -# Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de -# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des -# Solides) +# Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # -# Akantu 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. +# Akantu 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. # -# Akantu 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. +# Akantu 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 Akantu. If not, see . +# You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # #=============================================================================== add_mesh(test_grid_mesh circle.geo 2 1) register_test(test_grid test_grid.cc DEPENDS test_grid_mesh PACKAGE core) #register_test(test_types test_types.cc PACKAGE core) register_gtest_sources(SOURCES test_csr.cc PACKAGE core) register_gtest_sources(SOURCES test_arange_iterator.cc PACKAGE core HEADER_ONLY) register_gtest_sources(SOURCES test_zip_iterator.cc PACKAGE core HEADER_ONLY) register_gtest_sources(SOURCES test_array.cc PACKAGE core) register_gtest_sources(SOURCES test_tensors.cc PACKAGE core) register_gtest_test(test_common) diff --git a/test/test_common/test_arange_iterator.cc b/test/test_common/test_arange_iterator.cc index 62b3d13f2..eebc171da 100644 --- a/test/test_common/test_arange_iterator.cc +++ b/test/test_common/test_arange_iterator.cc @@ -1,69 +1,71 @@ /** * @file test_arange_iterator.cc * - * @author Nicolas Richart + * @author Nicolas Richart * - * @date creation Fri Aug 11 2017 + * @date creation: Fri Jun 18 2010 + * @date last modification: Sun Dec 03 2017 * - * @brief A Documented file. + * @brief A Documented file. * * @section LICENSE * - * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ + /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; TEST(TestArangeIterator, Stop) { size_t ref_i = 0; for (auto i : arange(10)) { EXPECT_EQ(ref_i, i); ++ref_i; } } TEST(TestArangeIterator, StartStop) { size_t ref_i = 1; for (auto i : arange(1, 10)) { EXPECT_EQ(ref_i, i); ++ref_i; } } TEST(TestArangeIterator, StartStopStep) { size_t ref_i = 1; for (auto i : arange(1, 22, 2)) { EXPECT_EQ(ref_i, i); ref_i += 2; } } TEST(TestArangeIterator, StartStopStepZipped) { int ref_i1 = -1, ref_i2 = 1; for (auto && i : zip(arange(-1, -10, -1), arange(1, 18, 2))) { EXPECT_EQ(ref_i1, std::get<0>(i)); EXPECT_EQ(ref_i2, std::get<1>(i)); ref_i1 += -1; ref_i2 += 2; } } diff --git a/test/test_common/test_array.cc b/test/test_common/test_array.cc index 60f1b20ee..6b587b845 100644 --- a/test/test_common/test_array.cc +++ b/test/test_common/test_array.cc @@ -1,255 +1,285 @@ +/** + * @file test_array.cc + * + * @author Nicolas Richart + * + * @date creation: Thu Nov 09 2017 + * @date last modification: Fri Jan 26 2018 + * + * @brief Test the arry class + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { class NonTrivial { public: NonTrivial() = default; NonTrivial(int a) : a(a){}; bool operator==(const NonTrivial & rhs) { return a == rhs.a; } int a{0}; }; bool operator==(const int & a, const NonTrivial & rhs) { return a == rhs.a; } std::ostream & operator<<(std::ostream & stream, const NonTrivial & _this) { stream << _this.a; return stream; } /* -------------------------------------------------------------------------- */ using TestTypes = ::testing::Types; /* -------------------------------------------------------------------------- */ ::testing::AssertionResult AssertType(const char * /*a_expr*/, const char * /*b_expr*/, const std::type_info & a, const std::type_info & b) { if (std::type_index(a) == std::type_index(b)) return ::testing::AssertionSuccess(); return ::testing::AssertionFailure() << debug::demangle(a.name()) << " != " << debug::demangle(b.name()) << ") are different"; } /* -------------------------------------------------------------------------- */ template class ArrayConstructor : public ::testing::Test { protected: using type = T; void SetUp() override { type_str = debug::demangle(typeid(T).name()); } template decltype(auto) construct(P &&... params) { return std::make_unique>(std::forward

(params)...); } protected: std::string type_str; }; TYPED_TEST_CASE(ArrayConstructor, TestTypes); TYPED_TEST(ArrayConstructor, ConstructDefault1) { auto array = this->construct(); EXPECT_EQ(0, array->size()); EXPECT_EQ(1, array->getNbComponent()); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault2) { auto array = this->construct(1000); EXPECT_EQ(1000, array->size()); EXPECT_EQ(1, array->getNbComponent()); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault3) { auto array = this->construct(1000, 10); EXPECT_EQ(1000, array->size()); EXPECT_EQ(10, array->getNbComponent()); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault4) { auto array = this->construct(1000, 10, "test"); EXPECT_EQ(1000, array->size()); EXPECT_EQ(10, array->getNbComponent()); EXPECT_STREQ("test", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault5) { auto array = this->construct(1000, 10, 1); EXPECT_EQ(1000, array->size()); EXPECT_EQ(10, array->getNbComponent()); EXPECT_EQ(1, array->operator()(10, 6)); EXPECT_STREQ("", array->getID().c_str()); } TYPED_TEST(ArrayConstructor, ConstructDefault6) { typename TestFixture::type defaultv[2] = {0, 1}; auto array = this->construct(1000, 2, defaultv); EXPECT_EQ(1000, array->size()); EXPECT_EQ(2, array->getNbComponent()); EXPECT_EQ(1, array->operator()(10, 1)); EXPECT_EQ(0, array->operator()(603, 0)); EXPECT_STREQ("", array->getID().c_str()); } /* -------------------------------------------------------------------------- */ template class ArrayFixture : public ArrayConstructor { public: void SetUp() override { ArrayConstructor::SetUp(); array = this->construct(1000, 10); } void TearDown() override { array.reset(nullptr); } protected: std::unique_ptr> array; }; TYPED_TEST_CASE(ArrayFixture, TestTypes); TYPED_TEST(ArrayFixture, Copy) { Array copy(*this->array); EXPECT_EQ(1000, copy.size()); EXPECT_EQ(10, copy.getNbComponent()); EXPECT_NE(this->array->storage(), copy.storage()); } TYPED_TEST(ArrayFixture, Set) { auto & arr = *(this->array); arr.set(12); EXPECT_EQ(12, arr(156, 5)); EXPECT_EQ(12, arr(520, 7)); EXPECT_EQ(12, arr(999, 9)); } TYPED_TEST(ArrayFixture, Resize) { auto & arr = *(this->array); auto ptr = arr.storage(); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_EQ(ptr, arr.storage()); EXPECT_EQ(1000, arr.getAllocatedSize()); arr.resize(3000); EXPECT_EQ(3000, arr.size()); EXPECT_EQ(3000, arr.getAllocatedSize()); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_EQ(nullptr, arr.storage()); EXPECT_EQ(0, arr.getAllocatedSize()); } TYPED_TEST(ArrayFixture, PushBack) { auto & arr = *(this->array); auto ptr = arr.storage(); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_EQ(ptr, arr.storage()); EXPECT_EQ(1000, arr.getAllocatedSize()); arr.resize(3000); EXPECT_EQ(3000, arr.size()); EXPECT_EQ(3000, arr.getAllocatedSize()); arr.resize(0); EXPECT_EQ(0, arr.size()); EXPECT_EQ(nullptr, arr.storage()); EXPECT_EQ(0, arr.getAllocatedSize()); } TYPED_TEST(ArrayFixture, ViewVector) { auto && view = make_view(*this->array, 10); EXPECT_NO_THROW(view.begin()); { auto it = view.begin(); EXPECT_EQ(10, it->size()); EXPECT_PRED_FORMAT2(AssertType, typeid(*it), typeid(Vector)); EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), typeid(VectorProxy)); } } TYPED_TEST(ArrayFixture, ViewMatrix) { { auto && view = make_view(*this->array, 2, 5); EXPECT_NO_THROW(view.begin()); { auto it = view.begin(); EXPECT_EQ(10, it->size()); EXPECT_EQ(2, it->size(0)); EXPECT_EQ(5, it->size(1)); EXPECT_PRED_FORMAT2(AssertType, typeid(*it), typeid(Matrix)); EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), typeid(MatrixProxy)); } } } TYPED_TEST(ArrayFixture, ViewVectorWrong) { auto && view = make_view(*this->array, 11); EXPECT_THROW(view.begin(), debug::ArrayException); } TYPED_TEST(ArrayFixture, ViewMatrixWrong) { auto && view = make_view(*this->array, 3, 7); EXPECT_THROW(view.begin(), debug::ArrayException); } TYPED_TEST(ArrayFixture, ViewMatrixIter) { std::size_t count = 0; for (auto && mat : make_view(*this->array, 10, 10)) { EXPECT_EQ(100, mat.size()); EXPECT_EQ(10, mat.size(0)); EXPECT_EQ(10, mat.size(1)); EXPECT_PRED_FORMAT2(AssertType, typeid(mat), typeid(Matrix)); ++count; } EXPECT_EQ(100, count); } TYPED_TEST(ArrayFixture, ConstViewVector) { const auto & carray = *this->array; auto && view = make_view(carray, 10); EXPECT_NO_THROW(view.begin()); { auto it = view.begin(); EXPECT_EQ(10, it->size()); EXPECT_PRED_FORMAT2(AssertType, typeid(*it), typeid(Vector)); EXPECT_PRED_FORMAT2(AssertType, typeid(it[0]), typeid(VectorProxy)); } } } // namespace diff --git a/test/test_common/test_csr.cc b/test/test_common/test_csr.cc index cf032e5fe..31f8aa9bc 100644 --- a/test/test_common/test_csr.cc +++ b/test/test_common/test_csr.cc @@ -1,103 +1,102 @@ /** * @file test_csr.cc * * @author Nicolas Richart * * @date creation: Mon Jul 30 2012 - * @date last modification: Sun Oct 19 2014 + * @date last modification: Sun Dec 03 2017 * * @brief Test the CSR (compressed sparse row) data structure * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_csr.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; class TestCsrFixture : public ::testing::Test { protected: void SetUp() override { csr.resizeRows(N); csr.clearRows(); for (UInt i = 0; i < N; ++i) { UInt nb_cols(UInt(rand() * double(N) / (RAND_MAX + 1.))); nb_cols_per_row.push_back(nb_cols); for (UInt j = 0; j < nb_cols; ++j) { ++csr.rowOffset(i); } } csr.countToCSR(); csr.resizeCols(); csr.beginInsertions(); for (UInt i = 0; i < N; ++i) { UInt nb_cols = nb_cols_per_row[i]; for (UInt j = 0; j < nb_cols; ++j) { csr.insertInRow(i, nb_cols - j); } } csr.endInsertions(); } std::vector nb_cols_per_row; CSR csr; size_t N = 1000; }; TEST_F(TestCsrFixture, CheckInsertion) { EXPECT_EQ(N, this->csr.getNbRows()); } TEST_F(TestCsrFixture, Iteration) { for (UInt i = 0; i < this->csr.getNbRows(); ++i) { auto it = this->csr.begin(i); auto end = this->csr.end(i); UInt nb_cols = this->nb_cols_per_row[i]; for (; it != end; ++it) { EXPECT_EQ(nb_cols, *it); nb_cols--; } EXPECT_EQ(0, nb_cols); } } TEST_F(TestCsrFixture, ReverseIteration) { for (UInt i = 0; i < csr.getNbRows(); ++i) { auto it = csr.rbegin(i); auto end = csr.rend(i); UInt nb_cols = nb_cols_per_row[i]; UInt j = nb_cols; for (; it != end; --it) { EXPECT_EQ((nb_cols - j + 1), *it); j--; } EXPECT_EQ(0, j); } } diff --git a/test/test_common/test_grid.cc b/test/test_common/test_grid.cc index d9093ca7f..b3ec583d3 100644 --- a/test/test_common/test_grid.cc +++ b/test/test_common/test_grid.cc @@ -1,108 +1,107 @@ /** * @file test_grid.cc * * @author Nicolas Richart * * @date creation: Thu Jul 15 2010 - * @date last modification: Thu Aug 06 2015 + * @date last modification: Fri Dec 08 2017 * * @brief Test the grid object * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_grid_dynamic.hh" #include "mesh.hh" #include "mesh_io.hh" using namespace akantu; int main(int argc, char * argv[]) { const UInt spatial_dimension = 2; akantu::initialize(argc, argv); Mesh circle(spatial_dimension); circle.read("circle.msh"); const auto & l = circle.getLocalLowerBounds(); const auto & u = circle.getLocalUpperBounds(); Real spacing[spatial_dimension] = {0.2, 0.2}; Vector s(spacing, spatial_dimension); Vector c = u; c += l; c /= 2.; SpatialGrid grid(spatial_dimension, s, c); Vector bary(spatial_dimension); Element el; el.ghost_type = _not_ghost; auto it = circle.firstType(spatial_dimension); auto last_type = circle.lastType(spatial_dimension); for (; it != last_type; ++it) { UInt nb_element = circle.getNbElement(*it); el.type = *it; for (UInt e = 0; e < nb_element; ++e) { el.element = e; circle.getBarycenter(el, bary); grid.insert(el, bary); } } std::cout << grid << std::endl; Mesh mesh(spatial_dimension, "save"); grid.saveAsMesh(mesh); mesh.write("grid.msh"); Vector pos(spatial_dimension); // const SpatialGrid::CellID & id = grid.getCellID(pos); // #if !defined AKANTU_NDEBUG // SpatialGrid::neighbor_cells_iterator nit = // grid.beginNeighborCells(id); // SpatialGrid::neighbor_cells_iterator nend = // grid.endNeighborCells(id); // for(;nit != nend; ++nit) { // std::cout << std::endl; // const SpatialGrid::Cell & cell = grid.getCell(*nit); // SpatialGrid::Cell::const_iterator cit = cell.begin(); // SpatialGrid::Cell::position_iterator pit = cell.begin_pos(); // SpatialGrid::Cell::const_iterator cend = cell.end(); // for (; cit != cend; ++cit, ++pit) { // std::cout << *cit << " " << *pit << std::endl; // } // } // #endif akantu::finalize(); return EXIT_SUCCESS; } diff --git a/test/test_common/test_tensors.cc b/test/test_common/test_tensors.cc index ec1797a10..f01a832a7 100644 --- a/test/test_common/test_tensors.cc +++ b/test/test_common/test_tensors.cc @@ -1,533 +1,563 @@ +/** + * @file test_tensors.cc + * + * @author Nicolas Richart + * + * @date creation: Tue Nov 14 2017 + * @date last modification: Mon Jan 22 2018 + * + * @brief test the tensors types + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { /* -------------------------------------------------------------------------- */ class TensorConstructorFixture : public ::testing::Test { public: void SetUp() override { for (auto & r : reference) { r = rand(); // google-test seeds srand() } } void TearDown() override {} template void compareToRef(const V & v) { for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i], v.storage()[i]); } } protected: const int size_{24}; const std::array mat_size{{4, 6}}; // const std::array tens3_size{{4, 2, 3}}; std::array reference; }; /* -------------------------------------------------------------------------- */ class TensorFixture : public TensorConstructorFixture { public: TensorFixture() : vref(reference.data(), size_), mref(reference.data(), mat_size[0], mat_size[1]) {} protected: Vector vref; Matrix mref; }; /* -------------------------------------------------------------------------- */ // Vector ---------------------------------------------------------------------- TEST_F(TensorConstructorFixture, VectorDefaultConstruct) { Vector v; EXPECT_EQ(0, v.size()); EXPECT_EQ(nullptr, v.storage()); EXPECT_EQ(false, v.isWrapped()); } TEST_F(TensorConstructorFixture, VectorConstruct1) { double r = rand(); Vector v(size_, r); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(r, v(i)); EXPECT_DOUBLE_EQ(r, v[i]); } } TEST_F(TensorConstructorFixture, VectorConstructWrapped) { Vector v(reference.data(), size_); EXPECT_EQ(size_, v.size()); EXPECT_EQ(true, v.isWrapped()); for (int i = 0; i < size_; ++i) { EXPECT_DOUBLE_EQ(reference[i], v(i)); EXPECT_DOUBLE_EQ(reference[i], v[i]); } } TEST_F(TensorConstructorFixture, VectorConstructInitializer) { Vector v{0., 1., 2., 3., 4., 5.}; EXPECT_EQ(6, v.size()); EXPECT_EQ(false, v.isWrapped()); for (int i = 0; i < 6; ++i) { EXPECT_DOUBLE_EQ(i, v(i)); } } TEST_F(TensorConstructorFixture, VectorConstructCopy1) { Vector vref(reference.data(), reference.size()); Vector v(vref); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructCopy2) { Vector vref(reference.data(), reference.size()); Vector v(vref, false); EXPECT_EQ(size_, v.size()); EXPECT_EQ(true, v.isWrapped()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructProxy1) { VectorProxy vref(reference.data(), reference.size()); EXPECT_EQ(size_, vref.size()); compareToRef(vref); Vector v(vref); EXPECT_EQ(size_, v.size()); EXPECT_EQ(true, v.isWrapped()); compareToRef(v); } TEST_F(TensorConstructorFixture, VectorConstructProxy2) { Vector vref(reference.data(), reference.size()); VectorProxy v(vref); EXPECT_EQ(size_, v.size()); compareToRef(v); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorEqual) { Vector v; v = vref; compareToRef(v); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); } TEST_F(TensorFixture, VectorEqualProxy) { VectorProxy vref_proxy(vref); Vector v; v = vref; compareToRef(v); EXPECT_EQ(size_, v.size()); EXPECT_EQ(false, v.isWrapped()); } TEST_F(TensorFixture, VectorEqualProxy2) { Vector v_store(size_, 0.); VectorProxy v(v_store); v = vref; compareToRef(v); compareToRef(v_store); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorSet) { Vector v(vref); compareToRef(v); double r = rand(); v.set(r); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(r, v[i]); } TEST_F(TensorFixture, VectorClear) { Vector v(vref); compareToRef(v); v.clear(); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0, v[i]); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, VectorDivide) { Vector v; double r = rand(); v = vref / r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, v[i]); } TEST_F(TensorFixture, VectorMultiply1) { Vector v; double r = rand(); v = vref * r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } TEST_F(TensorFixture, VectorMultiply2) { Vector v; double r = rand(); v = r * vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } TEST_F(TensorFixture, VectorAddition) { Vector v; v = vref + vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., v[i]); } TEST_F(TensorFixture, VectorSubstract) { Vector v; v = vref - vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., v[i]); } TEST_F(TensorFixture, VectorDivideEqual) { Vector v(vref); double r = rand(); v /= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, v[i]); } TEST_F(TensorFixture, VectorMultiplyEqual1) { Vector v(vref); double r = rand(); v *= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, v[i]); } TEST_F(TensorFixture, VectorMultiplyEqual2) { Vector v(vref); v *= v; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * reference[i], v[i]); } TEST_F(TensorFixture, VectorAdditionEqual) { Vector v(vref); v += vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., v[i]); } TEST_F(TensorFixture, VectorSubstractEqual) { Vector v(vref); v -= vref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., v[i]); } /* -------------------------------------------------------------------------- */ // Matrix ---------------------------------------------------------------------- TEST_F(TensorConstructorFixture, MatrixDefaultConstruct) { Matrix m; EXPECT_EQ(0, m.size()); EXPECT_EQ(0, m.rows()); EXPECT_EQ(0, m.cols()); EXPECT_EQ(nullptr, m.storage()); EXPECT_EQ(false, m.isWrapped()); } TEST_F(TensorConstructorFixture, MatrixConstruct1) { double r = rand(); Matrix m(mat_size[0], mat_size[1], r); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); for (int i = 0; i < mat_size[0]; ++i) { for (int j = 0; j < mat_size[1]; ++j) { EXPECT_EQ(r, m(i, j)); EXPECT_EQ(r, m[i + j * mat_size[0]]); } } } TEST_F(TensorConstructorFixture, MatrixConstructWrapped) { Matrix m(reference.data(), mat_size[0], mat_size[1]); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(true, m.isWrapped()); for (int i = 0; i < mat_size[0]; ++i) { for (int j = 0; j < mat_size[1]; ++j) { EXPECT_DOUBLE_EQ(reference[i + j * mat_size[0]], m(i, j)); } } compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructInitializer) { Matrix m{{0., 1., 2.}, {3., 4., 5.}}; EXPECT_EQ(6, m.size()); EXPECT_EQ(2, m.rows()); EXPECT_EQ(3, m.cols()); EXPECT_EQ(false, m.isWrapped()); int c = 0; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 3; ++j, ++c) { EXPECT_DOUBLE_EQ(c, m(i, j)); } } } TEST_F(TensorConstructorFixture, MatrixConstructCopy1) { Matrix mref(reference.data(), mat_size[0], mat_size[1]); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructCopy2) { Matrix mref(reference.data(), mat_size[0], mat_size[1]); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructProxy1) { MatrixProxy mref(reference.data(), mat_size[0], mat_size[1]); EXPECT_EQ(size_, mref.size()); EXPECT_EQ(mat_size[0], mref.size(0)); EXPECT_EQ(mat_size[1], mref.size(1)); compareToRef(mref); Matrix m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(true, m.isWrapped()); compareToRef(m); } TEST_F(TensorConstructorFixture, MatrixConstructProxy2) { Matrix mref(reference.data(), mat_size[0], mat_size[1]); MatrixProxy m(mref); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.size(0)); EXPECT_EQ(mat_size[1], m.size(1)); compareToRef(m); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixEqual) { Matrix m; m = mref; compareToRef(m); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); } TEST_F(TensorFixture, MatrixEqualProxy1) { MatrixProxy mref_proxy(mref); Matrix m; m = mref; compareToRef(m); EXPECT_EQ(size_, m.size()); EXPECT_EQ(mat_size[0], m.rows()); EXPECT_EQ(mat_size[1], m.cols()); EXPECT_EQ(false, m.isWrapped()); } TEST_F(TensorFixture, MatrixEqualProxy2) { Matrix m_store(mat_size[0], mat_size[1], 0.); MatrixProxy m(m_store); m = mref; compareToRef(m); compareToRef(m_store); } TEST_F(TensorFixture, MatrixEqualSlice) { Matrix m(mat_size[0], mat_size[1], 0.); for (unsigned int i = 0; i < m.cols(); ++i) m(i) = Vector(mref(i)); compareToRef(m); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixSet) { Matrix m(mref); compareToRef(m); double r = rand(); m.set(r); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(r, m[i]); } TEST_F(TensorFixture, MatrixClear) { Matrix m(mref); compareToRef(m); m.clear(); for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0, m[i]); } /* -------------------------------------------------------------------------- */ TEST_F(TensorFixture, MatrixDivide) { Matrix m; double r = rand(); m = mref / r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, m[i]); } TEST_F(TensorFixture, MatrixMultiply1) { Matrix m; double r = rand(); m = mref * r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, m[i]); } TEST_F(TensorFixture, MatrixMultiply2) { Matrix m; double r = rand(); m = r * mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, m[i]); } TEST_F(TensorFixture, MatrixAddition) { Matrix m; m = mref + mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., m[i]); } TEST_F(TensorFixture, MatrixSubstract) { Matrix m; m = mref - mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., m[i]); } TEST_F(TensorFixture, MatrixDivideEqual) { Matrix m(mref); double r = rand(); m /= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] / r, m[i]); } TEST_F(TensorFixture, MatrixMultiplyEqual1) { Matrix m(mref); double r = rand(); m *= r; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * r, m[i]); } TEST_F(TensorFixture, MatrixAdditionEqual) { Matrix m(mref); m += mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(reference[i] * 2., m[i]); } TEST_F(TensorFixture, MatrixSubstractEqual) { Matrix m(mref); m -= mref; for (int i = 0; i < size_; ++i) EXPECT_DOUBLE_EQ(0., m[i]); } #if defined(AKANTU_USE_LAPACK) TEST_F(TensorFixture, MatrixEigs) { Matrix m{{0, 1, 0, 0}, {1., 0, 0, 0}, {0, 1, 0, 1}, {0, 0, 4, 0}}; Matrix eig_vects(4, 4); Vector eigs(4); m.eig(eigs, eig_vects); Vector eigs_ref{2, 1., -1., -2}; auto lambda_v = m * eig_vects; for (int i = 0; i < 4; ++i) { EXPECT_NEAR(eigs_ref(i), eigs(i), 1e-14); for (int j = 0; j < 4; ++j) { EXPECT_NEAR(lambda_v(i)(j), eigs(i) * eig_vects(i)(j), 1e-14); } } } #endif /* -------------------------------------------------------------------------- */ } // namespace diff --git a/test/test_common/test_types.cc b/test/test_common/test_types.cc index 375768c12..38e180bf7 100644 --- a/test/test_common/test_types.cc +++ b/test/test_common/test_types.cc @@ -1,355 +1,355 @@ /** * @file test_types.cc * * @author Nicolas Richart * * @date creation: Fri May 15 2015 - * @date last modification: Sat Jul 11 2015 + * @date last modification: Wed Jun 14 2017 * * @brief Test the types declared in aka_types.hh * * @section LICENSE * - * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory - * (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include #include #include using namespace akantu; const Real tolerance = 1e-15; std::string itoa(UInt a) { std::stringstream sstr; sstr << a; return sstr.str(); } UInt testcounter = 0; struct wrap_error : std::runtime_error { wrap_error(const std::string & msg) : std::runtime_error(msg) {} }; struct size_error : std::runtime_error { size_error(const std::string & msg) : std::runtime_error(msg) {} }; struct data_error : std::runtime_error { data_error(const std::string & msg, UInt i) : std::runtime_error(msg), index(i) {} UInt index; }; template void compare_storages_with_ref(const type & a, Real * ref, UInt size, UInt line, const std::string & txt) { std::cout << std::setw(3) << (testcounter++) << ": " << std::setw(10) << txt << " - " << a << " - wrapped: " << std::boolalpha << a.isWrapped() << std::endl; if (a.size() != size) throw size_error("the size is not correct " + itoa(a.size()) + " instead of " + itoa(size) + " [Test at line: " + itoa(line) + "]"); Real * a_ptr = a.storage(); for (UInt i = 0; i < a.size(); ++i) { if (!((std::abs(a_ptr[i]) < tolerance && std::abs(ref[i]) < tolerance) || std::abs((a_ptr[i] - ref[i]) / a_ptr[i]) < tolerance)) { std::stringstream txt; txt << " std::abs(" << a_ptr[i] << " - " << ref[i] << " [= " << std::abs(a_ptr[i] - ref[i]) << "] ) > " << tolerance; throw data_error("storage differs at index " + itoa(i) + " [Test at line: " + itoa(line) + "]" + txt.str(), i); } } if (a_ptr == ref && !a.isWrapped()) throw wrap_error( "the storage should be wrapped but it is not [Test at line: " + itoa(line) + "]"); if (a_ptr != ref && a.isWrapped()) throw wrap_error( "the storage should not be wrapped but it is [Test at line: " + itoa(line) + "]"); } #define COMPARE(a, aref, txt) \ compare_storages_with_ref(a, aref, sizeof(aref) / sizeof(aref[0]), __LINE__, \ txt) #define COMPARE_STORAGE(a, aref, txt) \ compare_storages_with_ref(a, aref.storage(), aref.size(), __LINE__, txt) const UInt ref_size = 10; // clang-format off /* -------------------------------------------------------------------------- */ void test_constructor() { std::cout << "=== Test constructors ===" << std::endl; Real ref1[ref_size] = { 0. }; Real ref2[ref_size] = { 1563.58, 1563.58, 1563.58, 1563.58, 1563.58, 1563.58, 1563.58, 1563.58, 1563.58, 1563.58 }; Real ref3[ref_size] = { 23.1594, 79.6184, 77.9052, 47.9922, 12.8674, 37.1445, 64.8991, 80.3364, 98.4064, 73.7858 }; std::cout << "-- Vectors: " << std::endl; Vector v0 = { 23.1594, 79.6184, 77.9052, 47.9922, 12.8674, 37.1445, 64.8991, 80.3364, 98.4064, 73.7858 }; ; COMPARE ( v0, ref3, "init_list" ); Vector v1(ref_size); COMPARE ( v1, ref1, "normal" ); Vector v2(ref_size, 1563.58); COMPARE ( v2, ref2, "defval" ); Vector v3(ref3, ref_size); COMPARE ( v3, ref3, "wrapped" ); Vector v3dcw(v3); COMPARE ( v3dcw, ref3, "wdeepcopy" ); Vector v3scw(v3, false); COMPARE ( v3scw, ref3, "wshallow" ); Vector v3dc(v3dcw); COMPARE_STORAGE( v3dc, v3dcw, "deepcopy" ); Vector v3sc(v3dcw, false); COMPARE_STORAGE( v3sc, v3dcw, "shallow" ); VectorProxy vp1(ref3, ref_size); Vector v4(vp1); COMPARE ( v4, ref3, "proxyptr" ); VectorProxy vp2(v3dcw); Vector v5(vp2); COMPARE_STORAGE( v5, v3dcw, "proxyvdc" ); VectorProxy vp3(v3scw); Vector v6(vp3); COMPARE ( v6, ref3, "proxyvsc" ); /* ------------------------------------------------------------------------ */ std::cout << "-- Matrices: " << std::endl; Matrix m0 = {{23.1594, 37.1445}, {79.6184, 64.8991}, {77.9052, 80.3364}, {47.9922, 98.4064}, {12.8674, 73.7858}}; COMPARE ( m0, ref3 , "init_list" ); Matrix m1(5, 2); COMPARE ( m1, ref1 , "normal" ); Matrix m1t(2, 5); COMPARE ( m1t, ref1 , "tnormal" ); Matrix m2(5, 2, 1563.58); COMPARE ( m2, ref2 , "defval" ); Matrix m2t(2, 5, 1563.58); COMPARE ( m2t, ref2 , "tdefval" ); Matrix m3(ref3, 5, 2); COMPARE ( m3, ref3 , "wrapped" ); Matrix m3t(ref3, 2, 5); COMPARE ( m3t, ref3 , "twrapped" ); Matrix m3dcw(m3); COMPARE ( m3dcw, ref3 , "wdeepcopy" ); Matrix m3scw(m3, false); COMPARE ( m3scw, ref3 , "wshallow" ); Matrix m3dc(m3dcw); COMPARE_STORAGE( m3dc, m3dcw , "deepcopy" ); Matrix m3sc(m3dcw, false); COMPARE_STORAGE( m3sc, m3dcw , "shallow" ); Matrix m3tdcw(m3t); COMPARE (m3tdcw, ref3 , "twdeepcopy"); Matrix m3tscw(m3t, false); COMPARE (m3tscw, ref3 , "twshallow" ); Matrix m3tdc(m3tdcw); COMPARE_STORAGE( m3tdc, m3tdcw, "tdeepcopy" ); Matrix m3tsc(m3tdcw, false); COMPARE_STORAGE( m3tsc, m3tdcw, "tshallow" ); MatrixProxy mp1(ref3, 5, 2); Matrix m4(mp1); COMPARE ( m4, ref3, "proxyptr" ); MatrixProxy mp2(m3dcw); Matrix m5(mp2); COMPARE_STORAGE( m5, m3dcw, "proxyvdc" ); MatrixProxy mp3(m3scw); Matrix m6(mp3); COMPARE ( m6, ref3, "proxyvsc" ); MatrixProxy mp1t(ref3, 2, 5); Matrix m4t(mp1t); COMPARE ( m4t, ref3, "tproxyptr" ); MatrixProxy mp2t(m3tdcw); Matrix m5t(mp2t); COMPARE_STORAGE( m5t, m3tdcw, "tproxyvdc" ); MatrixProxy mp3t(m3tscw); Matrix m6t(mp3t); COMPARE ( m6t, ref3, "tproxyvsc" ); } /* -------------------------------------------------------------------------- */ void test_equal_and_accessors() { std::cout << "=== Test operator=() ===" << std::endl; Real ref[ref_size] = { 23.1594, 79.6184, 77.9052, 47.9922, 12.8674, 37.1445, 64.8991, 80.3364, 98.4064, 73.7858 }; Real mod[ref_size] = { 98.7982, 72.1227, 19.7815, 57.6722, 47.1088, 14.9865, 13.3171, 62.7973, 33.9493, 98.3052 }; std::cout << "-- Vectors: " << std::endl; Vector v (ref, ref_size); Vector vm(mod, ref_size); Vector vref1(v); Vector v1; v1 = vref1; COMPARE_STORAGE(v1, vref1, "simple=" ); for (UInt i = 0; i < ref_size; ++i) v1 (i) = mod[i]; COMPARE (v1, mod, "s_acces" ); COMPARE_STORAGE(vref1, v, "refcheck1"); Vector v2 = vref1; COMPARE_STORAGE(v2, vref1, "construc="); for (UInt i = 0; i < ref_size; ++i) v2 (i) = mod[i]; COMPARE (v2, mod, "c_acces" ); COMPARE_STORAGE(vref1, v, "refcheck2"); Vector vref2(vref1, false); Vector v1w; v1w = vref2; COMPARE_STORAGE(v1w, vref1, "w_simple=" ); for (UInt i = 0; i < ref_size; ++i) v1w(i) = mod[i]; COMPARE (v1w, mod, "ws_acces" ); try { COMPARE(vref2, ref, "refcheck3"); } catch(wrap_error &) {} Vector v2w = vref2; COMPARE_STORAGE(v2w, vref1, "w_constru="); for (UInt i = 0; i < ref_size; ++i) v2w(i) = mod[i]; COMPARE (v2w, mod, "wc_acces" ); try { COMPARE(vref2, ref, "refcheck4"); } catch(wrap_error &) {} VectorProxy vp1(vref1); Vector v3; v3 = vp1; COMPARE_STORAGE(v3, vref1, "p_simple=" ); for (UInt i = 0; i < ref_size; ++i) v3(i) = mod[i]; COMPARE (v3, mod, "ps_acces" ); COMPARE_STORAGE(vref1, v, "refcheck5"); Vector v4 = vp1; COMPARE_STORAGE(v4, vref1, "p_constru="); for (UInt i = 0; i < ref_size; ++i) v4(i) = mod[i]; try { COMPARE(v4, mod, "pc_acces" ); } catch (wrap_error &) {} COMPARE(vref1, mod, "refcheck6"); try { COMPARE(vref2, mod, "refcheck7"); } catch(wrap_error &) {} vref2 = v; VectorProxy vp2(vref2); Vector v3w; v3w = vp2; COMPARE_STORAGE(v3w, vref1, "pw_simpl="); for (UInt i = 0; i < ref_size; ++i) v3w(i) = mod[i]; COMPARE (v3w, mod, "pws_acces"); try { COMPARE(vref2, ref, "refcheck8"); } catch(wrap_error &) {} Vector v4w = vp2; COMPARE_STORAGE( v4w, vref1, "pw_constr="); for (UInt i = 0; i < ref_size; ++i) v4w(i) = mod[i]; try { COMPARE(v4w, mod, "pwc_acces"); } catch (wrap_error &) {} COMPARE_STORAGE(v4w, vref2, "refcheck9"); try { COMPARE(vref2, mod, "refcheck10"); } catch(wrap_error &) {} vref1 = v; Real store[ref_size] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; Vector vs(store, 10); VectorProxy vp3(vs); vp3 = vref1; try { COMPARE(vref1, store, "vp_equal_v"); } catch(wrap_error &) {} // Vector vref3(vm); // VectorProxy vp4 = vref3; // vp3 = vp4; // try { COMPARE(vs, mod, "vp_equal_vp"); } catch(wrap_error &) {} /* ------------------------------------------------------------------------ */ std::cout << "-- Matrices: " << std::endl; Matrix m (ref, 5, 2); Matrix mt(ref, 2, 5); Matrix m1 (5, 2); Matrix m1t(2, 5); for (UInt i = 0; i < 5; ++i) { for (UInt j = 0; j < 2; ++j) { m1(i, j) = ref[i + j*5]; m1t(j, i) = ref[j + i*2]; } } COMPARE_STORAGE( m1, m, "access" ); COMPARE_STORAGE(m1t, m, "t_access"); Matrix mm (mod, 5, 2); Matrix mmt(mod, 2, 5); Matrix m2(m); Matrix m3(m); for (UInt j = 0; j < 2; ++j) { Vector v = m2(j); for (UInt i = 0; i < 5; ++i) v(i) = mm(i, j); } COMPARE_STORAGE(m2, mm, "slicing"); for (UInt j = 0; j < 2; ++j) m3(j) = mm(j); COMPARE_STORAGE(m3, mm, "slic_slic"); COMPARE(mm, mod, "refcheck"); Real mod_1[ref_size] = { 98.7982, 72.1227, 197.815, 57.6722, 47.1088, 14.9865, 13.3171, 627.973, 33.9493, 98.3052 }; Matrix m4 (mm); m4 (2,0) = 197.815; m4 (2,1) = 627.973; COMPARE(m4, mod_1, "partial"); Matrix m4t(mmt); m4t(0,1) = 197.815; m4t(1,3) = 627.973; COMPARE(m4t, mod_1, "t_partial"); } /* -------------------------------------------------------------------------- */ void test_simple_operators() { std::cout << "=== Test simple operation ===" << std::endl; Real ref[ref_size] = { 23.1594, 79.6184, 77.9052, 47.9922, 12.8674, 37.1445, 64.8991, 80.3364, 98.4064, 73.7858 }; Real mod[ref_size] = { 98.7982, 72.1227, 19.7815, 57.6722, 47.1088, 14.9865, 13.3171, 62.7973, 33.9493, 98.3052 }; Real ref_div[ref_size] = { 1.163905920192984e+00, 4.001326766509196e+00, 3.915227661071464e+00, 2.411910744798472e+00, 6.466680068348578e-01, 1.866745401547894e+00, 3.261589104432606e+00, 4.037410795054780e+00, 4.945542265554328e+00, 3.708201829329581e+00 }; Real ref_tim[ref_size] = { 4.608257412000000e+02, 1.584246923200000e+03, 1.550157669600000e+03, 9.549487955999999e+02, 2.560355252000000e+02, 7.391012610000000e+02, 1.291362291800000e+03, 1.598533687200000e+03, 1.958090547200000e+03, 1.468189848400000e+03 }; Real ref_p_mod[ref_size] = { 1.219576000000000e+02, 1.517411000000000e+02, 9.768670000000000e+01, 1.056644000000000e+02, 5.997620000000001e+01, 5.213100000000000e+01, 7.821620000000000e+01, 1.431337000000000e+02, 1.323557000000000e+02, 1.720910000000000e+02 }; Real ref_m_mod[ref_size] = { -7.563879999999999e+01, 7.495699999999999e+00, 5.812369999999999e+01, -9.680000000000000e+00, -3.424140000000000e+01, 2.215800000000000e+01, 5.158200000000001e+01, 1.753910000000000e+01, 6.445710000000000e+01, -2.451940000000000e+01 }; std::cout << "-- Vectors: " << std::endl; Vector v (ref, ref_size); Vector vm(mod, ref_size); Vector vref(v); Vector vmod(vm); Vector v1 = vref / 19.898; COMPARE(v1, ref_div, "v / s" ); Vector v2 = vref * 19.898; COMPARE(v2, ref_tim, "v * s" ); Vector v3 = 19.898 * vref; COMPARE(v3, ref_tim, "s * v" ); Vector v4 = vref + vmod; COMPARE(v4, ref_p_mod, "v1 + v2" ); Vector v5 = vref - vmod; COMPARE(v5, ref_m_mod, "v1 - v2" ); Vector v6 = vref; v6 *= 19.898; COMPARE(v6, ref_tim, "v *= s" ); Vector v7 = vref; v7 /= 19.898; COMPARE(v7, ref_div, "v /= s" ); Vector v8 = vref; v8 += vmod; COMPARE(v8, ref_p_mod, "v1 += v2"); Vector v9 = vref; v9 -= vmod; COMPARE(v9, ref_m_mod, "v1 -= v2"); std::cout << "-- Matrices: " << std::endl; Matrix m (ref, 5, 2); Matrix mm(mod, 5, 2); Matrix mref(m); Matrix mmod(mm); Matrix m1 = mref / 19.898; COMPARE(m1, ref_div, "m / s" ); Matrix m2 = mref * 19.898; COMPARE(m2, ref_tim, "m * s" ); Matrix m3 = 19.898 * mref; COMPARE(m3, ref_tim, "s * m" ); Matrix m4 = mref + mmod; COMPARE(m4, ref_p_mod, "m1 + m2" ); Matrix m5 = mref - mmod; COMPARE(m5, ref_m_mod, "m1 - m2" ); Matrix m6 = mref; m6 *= 19.898; COMPARE(m6, ref_tim, "m *= s" ); Matrix m7 = mref; m7 /= 19.898; COMPARE(m7, ref_div, "m /= s" ); Matrix m8 = mref; m8 += mmod; COMPARE(m8, ref_p_mod, "m1 += m2"); Matrix m9 = mref; m9 -= mmod; COMPARE(m9, ref_m_mod, "m1 -= m2"); } // clang-format on /* -------------------------------------------------------------------------- */ int main() { test_constructor(); test_equal_and_accessors(); test_simple_operators(); return 0; } diff --git a/test/test_common/test_zip_iterator.cc b/test/test_common/test_zip_iterator.cc index d5957fb37..e8a39e0a7 100644 --- a/test/test_common/test_zip_iterator.cc +++ b/test/test_common/test_zip_iterator.cc @@ -1,179 +1,181 @@ /** * @file test_zip_iterator.cc * - * @author Nicolas Richart + * @author Nicolas Richart * - * @date creation Fri Jul 21 2017 + * @date creation: Fri Jul 21 2017 + * @date last modification: Fri Dec 08 2017 * - * @brief test the zip container and iterator + * @brief test the zip container and iterator * * @section LICENSE * - * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ + /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; template class A { public: A() = default; A(T a) : a(a){}; A(const A & other) : a(other.a), copy_counter(other.copy_counter + 1), move_counter(other.move_counter) {} A & operator=(const A & other) { if (this != &other) { a = other.a; copy_counter = other.copy_counter + 1; } return *this; } A(A && other) : a(std::move(other.a)), copy_counter(std::move(other.copy_counter)), move_counter(std::move(other.move_counter) + 1) {} A & operator=(A && other) { if (this != &other) { a = std::move(other.a); copy_counter = std::move(other.copy_counter); move_counter = std::move(other.move_counter) + 1; } return *this; } A & operator*=(const T & b) { a *= b; return *this; } T a; size_t copy_counter{0}; size_t move_counter{0}; }; template struct C { struct iterator { using reference = A; using difference_type = void; using iterator_category = std::input_iterator_tag; using value_type = A; using pointer = A *; iterator(T pos) : pos(std::move(pos)) {} A operator*() { return A(pos); } bool operator!=(const iterator & other) const { return pos != other.pos; } bool operator==(const iterator & other) const { return pos == other.pos; } iterator & operator++() { ++pos; return *this; } T pos; }; C(T begin_, T end_) : begin_(std::move(begin_)), end_(std::move(end_)) {} iterator begin() { return iterator(begin_); } iterator end() { return iterator(end_); } T begin_, end_; }; class TestZipFixutre : public ::testing::Test { protected: void SetUp() override { a.reserve(size); b.reserve(size); for (size_t i = 0; i < size; ++i) { a.emplace_back(i); b.emplace_back(i + size); } } template void check(A && a, B && b, size_t pos, size_t nb_copy, size_t nb_move) { EXPECT_EQ(pos, a.a); EXPECT_EQ(nb_copy, a.copy_counter); EXPECT_EQ(nb_move, a.move_counter); EXPECT_FLOAT_EQ(pos + this->size, b.a); EXPECT_EQ(nb_copy, b.copy_counter); EXPECT_EQ(nb_move, b.move_counter); } protected: size_t size{20}; std::vector> a; std::vector> b; }; TEST_F(TestZipFixutre, SimpleTest) { size_t i = 0; for (auto && pair : zip(this->a, this->b)) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 0); ++i; } } TEST_F(TestZipFixutre, ConstTest) { size_t i = 0; const auto & ca = this->a; const auto & cb = this->b; for (auto && pair : zip(ca, cb)) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 0); EXPECT_EQ(true, std::is_const< std::remove_reference_t(pair))>>::value); EXPECT_EQ(true, std::is_const< std::remove_reference_t(pair))>>::value); ++i; } } TEST_F(TestZipFixutre, MixteTest) { size_t i = 0; const auto & cb = this->b; for (auto && pair : zip(a, cb)) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 0); EXPECT_EQ(false, std::is_const< std::remove_reference_t(pair))>>::value); EXPECT_EQ(true, std::is_const< std::remove_reference_t(pair))>>::value); ++i; } } TEST_F(TestZipFixutre, MoveTest) { size_t i = 0; for (auto && pair : zip(C(0, this->size), C(this->size, 2 * this->size))) { this->check(std::get<0>(pair), std::get<1>(pair), i, 0, 1); ++i; } } diff --git a/test/test_fe_engine/CMakeLists.txt b/test/test_fe_engine/CMakeLists.txt index da791c77a..d2ca3a8b2 100644 --- a/test/test_fe_engine/CMakeLists.txt +++ b/test/test_fe_engine/CMakeLists.txt @@ -1,134 +1,127 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux +# @author Lucas Frerot +# @author Nicolas Richart # # @date creation: Fri Sep 03 2010 -# @date last modification: Mon Dec 07 2015 +# @date last modification: Fri Jan 26 2018 # # @brief configuration for FEM tests # # @section LICENSE # -# Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de -# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des -# Solides) +# Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # -# Akantu 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. +# Akantu 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. # -# Akantu 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. +# Akantu 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 Akantu. If not, see . +# You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== function(register_fem_test operation type) set(_target test_${operation}${type}) register_test(${_target} SOURCES test_${operation}.cc FILES_TO_COPY ${type}.msh COMPILE_OPTIONS TYPE=${type} PACKAGE core ) endfunction() #=============================================================================== macro(register_mesh_types package) package_get_element_types(${package} _types) foreach(_type ${_types}) if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${_type}.msh) list(APPEND _meshes ${_type}.msh) #register_fem_test(fe_engine_precomputation ${_type }) else() if(NOT ${_type} STREQUAL _point_1) message("The mesh ${_type}.msh is missing, the fe_engine test cannot be activated without it") endif() endif() endforeach() endmacro(register_mesh_types) set(_meshes) register_mesh_types(core) package_is_activated(structural_mechanics has_structural_mechanics) if(has_structural_mechanics) register_mesh_types(structural_mechanics) endif() #add_mesh(test_fem_circle_1_mesh circle.geo 2 1 OUTPUT circle1.msh) #add_mesh(test_fem_circle_2_mesh circle.geo 2 2 OUTPUT circle2.msh) # Tests for class MeshData macro(register_typed_test test_name type value1 value2) set(target test_${test_name}_${type}) register_test(${target} SOURCES test_${test_name}.cc COMPILE_OPTIONS "TYPE=${type};VALUE1=${value1};VALUE2=${value2}" PACKAGE core ) endmacro() register_typed_test(mesh_data string \"5\" \"10\") register_typed_test(mesh_data UInt 5 10) add_mesh(test_boundary_msh cube.geo 3 1) add_mesh(test_boundary_msh_physical_names cube_physical_names.geo 3 1) register_test(test_mesh_boundary SOURCES test_mesh_boundary.cc DEPENDS test_boundary_msh test_boundary_msh_physical_names PACKAGE core) register_test(test_facet_element_mapping SOURCES test_facet_element_mapping.cc DEPENDS test_boundary_msh_physical_names PACKAGE core) akantu_pybind11_add_module(aka_test MODULE pybind11_akantu.cc) register_gtest_sources( SOURCES test_fe_engine_precomputation.cc PACKAGE core pybind11 DEPENDS aka_test ) register_gtest_sources( SOURCES test_fe_engine_precomputation_structural.cc PACKAGE structural_mechanics ) register_gtest_sources( SOURCES test_fe_engine_gauss_integration.cc PACKAGE core ) register_gtest_sources( SOURCES test_gradient.cc PACKAGE core ) register_gtest_sources( SOURCES test_integrate.cc PACKAGE core ) register_gtest_sources( SOURCES test_inverse_map.cc PACKAGE core ) register_gtest_test(test_fe_engine FILES_TO_COPY ${_meshes}) diff --git a/test/test_fe_engine/py_engine/py_engine.py b/test/test_fe_engine/py_engine/py_engine.py index 943fd2353..37bd98ac5 100644 --- a/test/test_fe_engine/py_engine/py_engine.py +++ b/test/test_fe_engine/py_engine/py_engine.py @@ -1,344 +1,355 @@ #!/usr/bin/env python3 +# ------------------------------------------------------------------------------ +__author__ = "Nicolas Richart" +__copyright__ = "Copyright (C) 2016-2018, EPFL (Ecole Polytechnique Fédérale" \ + " de Lausanne) Laboratory (LSMS - Laboratoire de Simulation" \ + " en Mécanique des Solides)" +__credits__ = ["Nicolas Richart"] +__license__ = "L-GPLv3" +__maintainer__ = "Nicolas Richart" +__email__ = "nicolas.richart@epfl.ch" +# ------------------------------------------------------------------------------ + __all__ = ['Shapes'] import numpy as np import numpy.polynomial.polynomial as poly import aka_test class Shapes(object): NATURAL_COORDS = { (1, 'quadrangle'): np.array([[-1.], [1.], [0.]]), (2, 'quadrangle'): np.array([[-1., -1.], [ 1., -1.], [ 1., 1.], [-1., 1.], [ 0., -1.], [ 1., 0.], [ 0., 1.], [-1., 0.]]), (3, 'quadrangle'): np.array([[-1., -1., -1.], [ 1., -1., -1.], [ 1., 1., -1.], [-1., 1., -1.], [-1., -1., 1.], [ 1., -1., 1.], [ 1., 1., 1.], [-1., 1., 1.], [ 0., -1., -1.], [ 1., 0., -1.], [ 0., 1., -1.], [-1., 0., -1.], [-1., -1., 0.], [ 1., -1., 0.], [ 1., 1., 0.], [-1., 1., 0.], [ 0., -1., 1.], [ 1., 0., 1.], [ 0., 1., 1.], [-1., 0., 1.]]), (2, 'triangle'): np.array([[0., 0.], [1., 0.], [0., 1,], [.5, 0.], [.5, .5], [0., .5]]), (3, 'triangle'): np.array([[0., 0., 0.], [1., 0., 0.], [0., 1., 0.], [0., 0., 1.], [.5, 0., 0.], [.5, .5, 0.], [0., .5, 0.], [0., 0., .5], [.5, 0., .5], [0., .5, .5]]), (3, 'pentahedron'): np.array([[-1., 1., 0.], [-1., 0., 1.], [-1., 0., 0.], [ 1., 1., 0.], [ 1., 0., 1.], [ 1., 0., 0.], [-1., .5, .5], [-1., 0., .5], [-1., .5, 0.], [ 0., 1., 0.], [ 0., 0., 1.], [ 0., 0., 0.], [ 1., .5, .5], [ 1., 0., .5], [ 1., .5, 0.], [ 0., .5, .5], [ 0., 0., .5], [ 0., .5, 0.]]), } QUADRATURE_W = { (1, 'quadrangle', 1): np.array([2.]), (1, 'quadrangle', 2): np.array([1., 1.]), (2, 'triangle', 1): np.array([1./2.]), (2, 'triangle', 2): np.array([1., 1., 1.])/6., (3, 'triangle', 1): np.array([1./6.]), (3, 'triangle', 2): np.array([1., 1., 1., 1.])/24., (2, 'quadrangle', 1): np.array([1., 1., 1., 1.]), (2, 'quadrangle', 2): np.array([1., 1., 1., 1.]), (3, 'quadrangle', 1): np.array([1., 1., 1., 1., 1., 1., 1., 1.]), (3, 'quadrangle', 2): np.array([1., 1., 1., 1., 1., 1., 1., 1.]), (3, 'pentahedron', 1): np.array([1., 1., 1., 1., 1., 1.])/6., (3, 'pentahedron', 2): np.array([1., 1., 1., 1., 1., 1.])/6., } _tet_a = (5. - np.sqrt(5.))/20. _tet_b = (5. + 3.*np.sqrt(5.))/20. QUADRATURE_G = { (1, 'quadrangle', 1): np.array([[0.]]), (1, 'quadrangle', 2): np.array([[-1.], [1.]])/np.sqrt(3.), (2, 'triangle', 1): np.array([[1., 1.]])/3., (2, 'triangle', 2): np.array([[1./6., 1./6.], [2./3, 1./6], [1./6., 2./3.]]), (3, 'triangle', 1): np.array([[1., 1., 1.]])/4., (3, 'triangle', 2): np.array([[_tet_a, _tet_a, _tet_a], [_tet_b, _tet_a, _tet_a], [_tet_a, _tet_b, _tet_a], [_tet_a, _tet_a, _tet_b]]), (2, 'quadrangle', 1): np.array([[-1., -1.], [ 1., -1.], [-1., 1.], [ 1., 1.]])/np.sqrt(3.), (2, 'quadrangle', 2): np.array([[-1., -1.], [ 1., -1.], [-1., 1.], [ 1., 1.]])/np.sqrt(3.), (3, 'quadrangle', 1): np.array([[-1., -1., -1.], [ 1., -1., -1.], [-1., 1., -1.], [ 1., 1., -1.], [-1., -1., 1.], [ 1., -1., 1.], [-1., 1., 1.], [ 1., 1., 1.]])/np.sqrt(3.), (3, 'quadrangle', 2): np.array([[-1., -1., -1.], [ 1., -1., -1.], [-1., 1., -1.], [ 1., 1., -1.], [-1., -1., 1.], [ 1., -1., 1.], [-1., 1., 1.], [ 1., 1., 1.]])/np.sqrt(3.), # (3, 'pentahedron', 1): np.array([[-1./np.sqrt(3.), 0.5, 0.5], # [-1./np.sqrt(3.), 0. , 0.5], # [-1./np.sqrt(3.), 0.5, 0. ], # [ 1./np.sqrt(3.), 0.5, 0.5], # [ 1./np.sqrt(3.), 0. , 0.5], # [ 1./np.sqrt(3.), 0.5 ,0. ]]), (3, 'pentahedron', 1): np.array([[-1./np.sqrt(3.), 1./6., 1./6.], [-1./np.sqrt(3.), 2./3., 1./6.], [-1./np.sqrt(3.), 1./6., 2./3.], [ 1./np.sqrt(3.), 1./6., 1./6.], [ 1./np.sqrt(3.), 2./3., 1./6.], [ 1./np.sqrt(3.), 1./6., 2./3.]]), (3, 'pentahedron', 2): np.array([[-1./np.sqrt(3.), 1./6., 1./6.], [-1./np.sqrt(3.), 2./3., 1./6.], [-1./np.sqrt(3.), 1./6., 2./3.], [ 1./np.sqrt(3.), 1./6., 1./6.], [ 1./np.sqrt(3.), 2./3., 1./6.], [ 1./np.sqrt(3.), 1./6., 2./3.]]), } ELEMENT_TYPES = { '_segment_2': ('quadrangle', 1, 'lagrange', 1, 2), '_segment_3': ('quadrangle', 2, 'lagrange', 1, 3), '_triangle_3': ('triangle', 1, 'lagrange', 2, 3), '_triangle_6': ('triangle', 2, 'lagrange', 2, 6), '_quadrangle_4': ('quadrangle', 1, 'serendip', 2, 4), '_quadrangle_8': ('quadrangle', 2, 'serendip', 2, 8), '_tetrahedron_4': ('triangle', 1, 'lagrange', 3, 4), '_tetrahedron_10': ('triangle', 2, 'lagrange', 3, 10), '_pentahedron_6': ('pentahedron', 1, 'lagrange', 3, 6), '_pentahedron_15': ('pentahedron', 2, 'lagrange', 3, 15), '_hexahedron_8': ('quadrangle', 1, 'serendip', 3, 8), '_hexahedron_20': ('quadrangle', 2, 'serendip', 3, 20), } MONOMES = {(1, 'quadrangle'): np.array([[0], [1], [2], [3], [4], [5]]), (2, 'triangle'): np.array([[0, 0], # 1 [1, 0], [0, 1], # x y [2, 0], [1, 1], [0, 2]]), # x^2 x.y y^2 (2, 'quadrangle'): np.array([[0, 0], [1, 0], [1, 1], [0, 1], [2, 0], [2, 1], [1, 2], [0, 2]]), (3, 'triangle'): np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [2, 0, 0], [1, 1, 0], [0, 2, 0], [0, 1, 1], [0, 0, 2], [1, 0, 1]]), (3, 'quadrangle'): np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0], [1, 0, 1], [0, 1, 1], [1, 1, 1], [2, 0, 0], [0, 2, 0], [0, 0, 2], [2, 1, 0], [2, 0, 1], [2, 1, 1], [1, 2, 0], [0, 2, 1], [1, 2, 1], [1, 0, 2], [0, 1, 2], [1, 1, 2]]), } SHAPES = { (3, 'pentahedron', 1): np.array([ [[[ 0., 0.], [ 1., 0.]], [[ 0., 0.], [-1., 0.]]], [[[ 0., 1.], [ 0., 0.]], [[ 0., -1.], [ 0., 0.]]], [[[ 1., -1.], [-1., 0.]], [[-1., 1.], [ 1., 0.]]], [[[ 0., 0.], [ 1., 0.]], [[ 0., 0.], [ 1., 0.]]], [[[ 0., 1.], [ 0., 0.]], [[ 0., 1.], [ 0., 0.]]], [[[ 1., -1.], [-1., 0.]], [[ 1., -1.], [-1., 0.]]] ])/2., (3, 'pentahedron', 2): np.array([ # 0 [[[ 0. , 0. , 0. ], [-1. , 0. , 0. ], [ 1. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0.5, 0. , 0. ], [-1. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0.5, 0. , 0. ], [ 0. , 0. , 0. ]]], # 1 [[[ 0. , -1. , 1. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0.5, -1. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0.5, 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 2 [[[ 0. , -1. , 1. ], [-1. , 2. , 0. ], [ 1. , 0. , 0. ]], [[-0.5, 1.5, -1. ], [ 1.5, -2. , 0. ], [-1. , 0. , 0. ]], [[ 0.5, -0.5, 0. ], [-0.5, 0. , 0. ], [ 0. , 0. , 0. ]]], # 3 [[[ 0. , 0. , 0. ], [-1. , 0. , 0. ], [ 1. , 0. , 0. ]], [[ 0. , 0. , 0. ], [-0.5, 0. , 0. ], [ 1. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0.5, 0. , 0. ], [ 0. , 0. , 0. ]]], # 4 [[[ 0. , -1. , 1. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , -0.5, 1. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0.5, 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 5 [[[ 0. , -1. , 1. ], [-1. , 2. , 0. ], [ 1. , 0. , 0. ]], [[ 0.5, -1.5, 1. ], [-1.5, 2. , 0. ], [ 1. , 0. , 0. ]], [[ 0.5, -0.5, 0. ], [-0.5, 0. , 0. ], [ 0. , 0. , 0. ]]], # 6 [[[ 0. , 0. , 0. ], [ 0. , 2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , -2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 7 [[[ 0. , 2. , -2. ], [ 0. , -2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , -2. , 2. ], [ 0. , 2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 8 [[[ 0. , 0. , 0. ], [ 2. , -2. , 0. ], [-2. , 0. , 0. ]], [[ 0. , 0. , 0. ], [-2. , 2. , 0. ], [ 2. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 9 [[[ 0. , 0. , 0. ], [ 1. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [-1. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 10 [[[ 0. , 1. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , -1. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 11 [[[ 1. , -1. , 0. ], [-1. , 0. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]], [[-1. , 1. , 0. ], [ 1. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 12 [[[ 0. , 0. , 0. ], [ 0. , 2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 13 [[[ 0. , 2. , -2. ], [ 0. , -2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 2. , -2. ], [ 0. , -2. , 0. ], [ 0. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], # 14 [[[ 0. , 0. , 0. ], [ 2. , -2. , 0. ], [-2. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 2. , -2. , 0. ], [-2. , 0. , 0. ]], [[ 0. , 0. , 0. ], [ 0. , 0. , 0. ], [ 0. , 0. , 0. ]]], ])} def __init__(self, element): self._shape, self._order, self._inter_poly, self._dim, self._nnodes = self.ELEMENT_TYPES[element] self._ksi = self.NATURAL_COORDS[(self._dim, self._shape)][:self._nnodes] self._g = self.QUADRATURE_G[(self._dim, self._shape, self._order)] self._w = self.QUADRATURE_W[(self._dim, self._shape, self._order)] def polyval(self, x, p): if 1 == self._dim: return poly.polyval(x[0], p) if 2 == self._dim: return poly.polyval2d(x[0], x[1], p) if 3 == self._dim: return poly.polyval3d(x[0], x[1], x[2], p) def shape_from_monomes(self): momo = self.MONOMES[(self._dim, self._shape)][:self._nnodes] _shape = list(momo[0]) for s in range(len(_shape)): _shape[s] = max(momo[:,s])+1 self._poly_shape = tuple(_shape) self._monome = [] for m in momo: p = np.zeros(self._poly_shape) p[tuple(m)] = 1 self._monome.append(p) # evaluate polynomial constant for shapes _x = self._ksi _xe = np.zeros((self._nnodes, self._nnodes)) for n in range(self._nnodes): _xe[:,n] = [self.polyval(_x[n], m) for m in self._monome] _a = np.linalg.inv(_xe) _n = np.zeros((self._nnodes,) + self._monome[0].shape) # set shapes polynomials for n in range(self._nnodes): for m in range(len(self._monome)): _n[n] += _a[n, m] * self._monome[m] return _n def compute_shapes(self): if (self._dim, self._shape) in self.MONOMES: return self.shape_from_monomes() else: _n = self.SHAPES[(self._dim, self._shape, self._order)] self._poly_shape = _n[0].shape return _n def precompute(self, **kwargs): X = np.array(kwargs["X"], copy=False) nb_element = X.shape[0] X = X.reshape(nb_element, self._nnodes, self._dim) _x = self._ksi _n = self.compute_shapes() # sanity check on shapes for n in range(self._nnodes): for m in range(self._nnodes): v = self.polyval(_x[n], _n[m]) ve = 1. if n == m else 0. test = np.isclose(v, ve) if not test: raise Exception("Most probably an error in the shapes evaluation") # compute shapes derivatives _b = np.zeros((self._dim, self._nnodes,) + self._poly_shape) for d in range(self._dim): for n in range(self._nnodes): _der = poly.polyder(_n[n], axis=d) _mshape = np.array(self._poly_shape) _mshape[d] = _mshape[d] - _der.shape[d] _mshape = tuple(_mshape) _comp = np.zeros(_mshape) if 1 == self._dim: _bt = np.hstack((_der, _comp)) else: if 0 == d: _bt = np.vstack((_der, _comp)) if 1 == d: _bt = np.hstack((_der, _comp)) if 2 == d: _bt = np.dstack((_der, _comp)) _b[d, n] = _bt _nb_quads = len(self._g) _nq = np.zeros((_nb_quads, self._nnodes)) _bq = np.zeros((_nb_quads, self._dim, self._nnodes)) # evaluate shapes and shapes derivatives on gauss points for q in range(_nb_quads): _g = self._g[q] for n in range(self._nnodes): _nq[q, n] = self.polyval(_g, _n[n]) for d in range(self._dim): _bq[q, d, n] = self.polyval(_g, _b[d, n]) _j = np.array(kwargs['j'], copy=False).reshape((nb_element, _nb_quads)) _B = np.array(kwargs['B'], copy=False).reshape((nb_element, _nb_quads, self._nnodes, self._dim)) _N = np.array(kwargs['N'], copy=False).reshape((nb_element, _nb_quads, self._nnodes)) _Q = np.array(kwargs['Q'], copy=False) if np.linalg.norm(_Q - self._g.T) > 1e-15: raise Exception('Not using the same quadrature points') for e in range(nb_element): for q in range(_nb_quads): _J = np.matmul(_bq[q], X[e]) if(np.linalg.norm(_N[e, q] - _nq[q]) > 1e-10): print(f"{e},{q}") print(_N[e, q]) print(_nq[q]) _N[e, q] = _nq[q] _tmp = np.matmul(np.linalg.inv(_J), _bq[q]) _B[e, q] = _tmp.T _j[e, q] = np.linalg.det(_J) * self._w[q] diff --git a/test/test_fe_engine/pybind11_akantu.cc b/test/test_fe_engine/pybind11_akantu.cc index fd8e247a8..7e2d16fb2 100644 --- a/test/test_fe_engine/pybind11_akantu.cc +++ b/test/test_fe_engine/pybind11_akantu.cc @@ -1,73 +1,103 @@ +/** + * @file pybind11_akantu.cc + * + * @author Nicolas Richart + * + * @date creation: Fri Dec 22 2017 + * @date last modification: Thu Dec 28 2017 + * + * @brief tool to wrap akantu classes in python + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + /* -------------------------------------------------------------------------- */ #include "pybind11_akantu.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; namespace akantu { PYBIND11_MODULE(aka_test, m) { m.doc() = "Module for the tests of akantu"; py::class_>(m, "ArrayProxy", py::buffer_protocol()) .def_buffer([](ArrayProxy & a) { return py::buffer_info( a.storage(), /* Pointer to buffer */ sizeof(Real), /* Size of one scalar */ py::format_descriptor::format(), /* Python struct-style format descriptor */ 2, /* Number of dimensions */ {a.size(), a.getNbComponent()}, /* Buffer dimensions */ {sizeof(Real) * a.getNbComponent(), /* Strides (in bytes) for each index */ sizeof(Real)}); }) .def(py::init([](py::array & n) { /* Request a buffer descriptor from Python */ py::buffer_info info = n.request(); /* Some sanity checks ... */ if (info.format != py::format_descriptor::format()) throw std::runtime_error( "Incompatible format: expected a double array!"); if (info.ndim != 2) throw std::runtime_error("Incompatible buffer dimension!"); return std::make_unique>(static_cast(info.ptr), info.shape[0], info.shape[1]); })); py::class_>(m, "Matrix", py::buffer_protocol()) .def_buffer([](MatrixProxy & a) { return py::buffer_info( a.storage(), /* Pointer to buffer */ sizeof(Real), /* Size of one scalar */ py::format_descriptor::format(), /* Python struct-style format descriptor */ 2, /* Number of dimensions */ {a.size(0), a.size(1)}, /* Buffer dimensions */ {sizeof(Real), a.size(0) * sizeof(Real)} /* Strides (in bytes) for each index */ ); }) .def(py::init([](py::array & n) { /* Request a buffer descriptor from Python */ py::buffer_info info = n.request(); /* Some sanity checks ... */ if (info.format != py::format_descriptor::format()) throw std::runtime_error( "Incompatible format: expected a double array!"); if (info.ndim != 2) throw std::runtime_error("Incompatible buffer dimension!"); return std::make_unique>( static_cast(info.ptr), info.shape[0], info.shape[1]); })); } // Module aka test } // namespace akantu diff --git a/test/test_fe_engine/pybind11_akantu.hh b/test/test_fe_engine/pybind11_akantu.hh index 3602c1bfe..f5e801144 100644 --- a/test/test_fe_engine/pybind11_akantu.hh +++ b/test/test_fe_engine/pybind11_akantu.hh @@ -1,110 +1,140 @@ +/** + * @file pybind11_akantu.hh + * + * @author Nicolas Richart + * + * @date creation: Fri Dec 22 2017 + * @date last modification: Thu Dec 28 2017 + * + * @brief tool to wrap akantu classes in python + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + #include "aka_array.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PYBIND11_AKANTU_HH__ #define __AKANTU_PYBIND11_AKANTU_HH__ namespace akantu { template class ArrayProxy : public Array { public: ArrayProxy(T * wrapped_memory, UInt size = 0, UInt nb_component = 1, const ID & id = "") : Array(0, nb_component, id) { this->values = wrapped_memory; this->size_ = size; } ArrayProxy(const ArrayProxy & array) : Array(0, array.nb_component, array.id) { this->values = array.values; this->size_ = array.size_; } ArrayProxy(ArrayProxy && array) { this->nb_component = std::move(array.nb_component); this->values = std::move(array.values); this->size_ = std::move(array.size_); this->id = std::move(array.id); } ArrayProxy(Array & array) : Array(0, array.getNbComponent(), array.getID()) { this->values = array.storage(); this->size_ = array.size(); } ~ArrayProxy() { this->values = nullptr; this->size_ = 0; } void setNbComponent(UInt nb_component) { UInt new_size = this->size_ / nb_component; AKANTU_DEBUG_ASSERT( nb_component * new_size == this->nb_component * this->size_, nb_component << " is not valid as a new number of component for this array"); this->nb_component = nb_component; this->size_ = new_size; } void resize(UInt new_size) { AKANTU_DEBUG_ASSERT(this->size_ == new_size, "cannot resize a temporary vector"); } }; template decltype(auto) make_proxy(Array & array) { return ArrayProxy(array); } template decltype(auto) make_proxy(const Matrix & array) { return MatrixProxy(array); } } // namespace akantu // namespace pybind11 { // namespace detail { // template <> struct type_caster> { // public: // PYBIND11_TYPE_CASTER(akantu::ArrayProxy, _("ArrayProxy")); // bool load(handle, bool) { // // /* Extract PyObject from handle */ // // PyObject *source = src.ptr(); // // /* Try converting into a Python integer value */ // // PyObject *tmp = PyNumber_Long(source); // // if (!tmp) // // return false; // // /* Now try to convert into a C++ int */ // // value.long_value = PyLong_AsLong(tmp); // // Py_DECREF(tmp); // // /* Ensure return code was OK (to avoid out-of-range errors etc) */ // // return !(value.long_value == -1 && !PyErr_Occurred()); // return false; // } // static handle cast(akantu::ArrayProxy & src, // return_value_policy /* policy */, handle parent) { // constexpr ssize_t elem_size = sizeof(akantu::Real); // ssize_t nb_comp = src.getNbComponent(); // ssize_t size = src.size(); // std::cout << "(src.storage()) // << ">\n"; // auto a = array({size, nb_comp}, {elem_size * nb_comp, elem_size}, // src.storage(), handle()); // return a.release(); // } // }; // } // namespace detail // } // namespace pybind11 #endif /* __AKANTU_PYBIND11_AKANTU_HH__ */ diff --git a/test/test_fe_engine/test_facet_element_mapping.cc b/test/test_fe_engine/test_facet_element_mapping.cc index 4a0e06fc8..541853e8f 100644 --- a/test/test_fe_engine/test_facet_element_mapping.cc +++ b/test/test_fe_engine/test_facet_element_mapping.cc @@ -1,135 +1,135 @@ /** * @file test_facet_element_mapping.cc * * @author Dana Christen * * @date creation: Fri May 03 2013 - * @date last modification: Mon Jul 13 2015 + * @date last modification: Tue Nov 07 2017 * * @brief Test of the MeshData class * * @section LICENSE * - * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_error.hh" #include "mesh.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; using namespace std; int main(int argc, char * argv[]) { // Testing the subelement-to-element mappings UInt spatial_dimension(3); akantu::initialize(argc, argv); Mesh mesh(spatial_dimension, "my_mesh"); mesh.read("./cube_physical_names.msh"); typedef Array> ElemToSubelemMapping; typedef Array SubelemToElemMapping; std::cout << "ELEMENT-SUBELEMENT MAPPING:" << std::endl; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(spatial_dimension, *git); Mesh::type_iterator tend = mesh.lastType(spatial_dimension, *git); std::cout << " " << "Ghost type: " << *git << std::endl; for (; tit != tend; ++tit) { const SubelemToElemMapping & subelement_to_element = mesh.getSubelementToElement(*tit, *git); std::cout << " " << " " << "Element type: " << *tit << std::endl; std::cout << " " << " " << " " << "subelement_to_element:" << std::endl; subelement_to_element.printself(std::cout, 8); for (UInt i(0); i < subelement_to_element.size(); ++i) { std::cout << " "; for (UInt j(0); j < mesh.getNbFacetsPerElement(*tit); ++j) { if (subelement_to_element(i, j) != ElementNull) { std::cout << subelement_to_element(i, j); std::cout << ", "; } else { std::cout << "ElementNull" << ", "; } } std::cout << "for element " << i << std::endl; } } tit = mesh.firstType(spatial_dimension - 1, *git); tend = mesh.lastType(spatial_dimension - 1, *git); for (; tit != tend; ++tit) { const ElemToSubelemMapping & element_to_subelement = mesh.getElementToSubelement(*tit, *git); std::cout << " " << " " << "Element type: " << *tit << std::endl; std::cout << " " << " " << " " << "element_to_subelement:" << std::endl; element_to_subelement.printself(std::cout, 8); for (UInt i(0); i < element_to_subelement.size(); ++i) { const std::vector & vec = element_to_subelement(i); std::cout << " "; std::cout << "item " << i << ": [ "; if (vec.size() > 0) { for (UInt j(0); j < vec.size(); ++j) { if (vec[j] != ElementNull) { std::cout << vec[j] << ", "; } else { std::cout << "ElementNull" << ", "; } } } else { std::cout << "empty, "; } std::cout << "]" << ", " << std::endl; } std::cout << std::endl; } } return 0; } diff --git a/test/test_fe_engine/test_fe_engine_fixture.hh b/test/test_fe_engine/test_fe_engine_fixture.hh index 8e2ba7b74..43d00d230 100644 --- a/test/test_fe_engine/test_fe_engine_fixture.hh +++ b/test/test_fe_engine/test_fe_engine_fixture.hh @@ -1,83 +1,113 @@ +/** + * @file test_fe_engine_fixture.hh + * + * @author Nicolas Richart + * + * @date creation: Tue Nov 14 2017 + * @date last modification: Mon Feb 19 2018 + * + * @brief Fixture for feengine tests + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__ #define __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__ using namespace akantu; /// Generic class for FEEngine tests template class shape_t, ElementKind kind = _ek_regular> class TestFEMBaseFixture : public ::testing::Test { public: static constexpr const ElementType type = type_::value; static constexpr const size_t dim = ElementClass::getSpatialDimension(); using FEM = FEEngineTemplate; /// Setup reads mesh corresponding to element type and initializes an FEEngine void SetUp() override { const auto dim = this->dim; const auto type = this->type; mesh = std::make_unique(dim); std::stringstream meshfilename; meshfilename << type << ".msh"; this->readMesh(meshfilename.str()); lower = mesh->getLowerBounds(); upper = mesh->getUpperBounds(); nb_element = this->mesh->getNbElement(type); fem = std::make_unique(*mesh, dim, "my_fem"); nb_quadrature_points_total = GaussIntegrationElement::getNbQuadraturePoints() * nb_element; SCOPED_TRACE(aka::to_string(type)); } void TearDown() override { fem.reset(nullptr); mesh.reset(nullptr); } /// Should be reimplemented if further treatment of the mesh is needed virtual void readMesh(std::string file_name) { mesh->read(file_name); } protected: std::unique_ptr fem; std::unique_ptr mesh; UInt nb_element; UInt nb_quadrature_points_total; Vector lower; Vector upper; }; template class shape_t, ElementKind kind> constexpr const ElementType TestFEMBaseFixture::type; template class shape_t, ElementKind kind> constexpr const size_t TestFEMBaseFixture::dim; /* -------------------------------------------------------------------------- */ /// Base class for test with Lagrange FEEngine and regular elements template using TestFEMFixture = TestFEMBaseFixture; /* -------------------------------------------------------------------------- */ using types = gtest_list_t; TYPED_TEST_CASE(TestFEMFixture, types); #endif /* __AKANTU_TEST_FE_ENGINE_FIXTURE_HH__ */ diff --git a/test/test_fe_engine/test_fe_engine_gauss_integration.cc b/test/test_fe_engine/test_fe_engine_gauss_integration.cc index 068899371..7cf2edd01 100644 --- a/test/test_fe_engine/test_fe_engine_gauss_integration.cc +++ b/test/test_fe_engine/test_fe_engine_gauss_integration.cc @@ -1,155 +1,154 @@ /** - * @file test_fe_engine_precomputation.cc + * @file test_fe_engine_gauss_integration.cc * * @author Nicolas Richart * - * @date creation: Mon Jun 14 2010 - * @date last modification: Mon Jul 13 2015 + * @date creation: Tue May 24 2016 + * @date last modification: Mon Feb 19 2018 * - * @brief test integration on elements, this test consider that mesh is a cube + * @brief test integration on elements, this test consider that mesh is a cube * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { /* -------------------------------------------------------------------------- */ template using degree_t = std::integral_constant; /* -------------------------------------------------------------------------- */ using TestDegreeTypes = std::tuple, degree_t<1>, degree_t<2>, degree_t<3>, degree_t<4>, degree_t<5>>; std::array, 3> global_polys{ {{0.40062394, 0.13703225, 0.51731446, 0.87830084, 0.5410543, 0.71842292}, {0.41861835, 0.11080576, 0.49874043, 0.49077504, 0.85073835, 0.66259755}, {0.92620845, 0.7503478, 0.62962232, 0.31662719, 0.64069644, 0.30878135}}}; template class TestGaussIntegrationFixture : public TestFEMFixture> { protected: using parent = TestFEMFixture>; static constexpr size_t degree{std::tuple_element_t<1, T>::value}; public: TestGaussIntegrationFixture() : integration_points_pos(0, parent::dim) {} void SetUp() override { parent::SetUp(); this->fem->initShapeFunctions(); auto integration_points = this->fem->getIntegrator().template getIntegrationPoints < parent::type, degree == 0 ? 1 : degree > (); nb_integration_points = integration_points.cols(); auto shapes_size = ElementClass::getShapeSize(); Array shapes(0, shapes_size); this->fem->getShapeFunctions() .template computeShapesOnIntegrationPoints( this->mesh->getNodes(), integration_points, shapes, _not_ghost); auto vect_size = this->nb_integration_points * this->nb_element; integration_points_pos.resize(vect_size); this->fem->getShapeFunctions() .template interpolateOnIntegrationPoints( this->mesh->getNodes(), integration_points_pos, this->dim, shapes); for (size_t d = 0; d < this->dim; ++d) { polys[d] = global_polys[d].extract(degree); } } void testIntegrate() { std::stringstream sstr; sstr << this->type << ":" << this->degree; SCOPED_TRACE(sstr.str().c_str()); auto vect_size = this->nb_integration_points * this->nb_element; Array polynomial(vect_size); size_t dim = parent::dim; for (size_t d = 0; d < dim; ++d) { auto poly = this->polys[d]; for (auto && pair : zip(polynomial, make_view(this->integration_points_pos, dim))) { auto && p = std::get<0>(pair); auto & x = std::get<1>(pair); p = poly(x(d)); } auto res = this->fem->getIntegrator() .template integrate( polynomial); auto expect = poly.integrate(this->lower(d), this->upper(d)); for (size_t o = 0; o < dim; ++o) { if (o == d) continue; expect *= this->upper(d) - this->lower(d); } EXPECT_NEAR(expect, res, 5e-14); } } protected: UInt nb_integration_points; std::array, parent::dim> polynomial; Array integration_points_pos; std::array, 3> polys; }; template constexpr size_t TestGaussIntegrationFixture::degree; /* -------------------------------------------------------------------------- */ /* Tests */ /* -------------------------------------------------------------------------- */ TYPED_TEST_CASE_P(TestGaussIntegrationFixture); TYPED_TEST_P(TestGaussIntegrationFixture, ArbitraryOrder) { this->testIntegrate(); } REGISTER_TYPED_TEST_CASE_P(TestGaussIntegrationFixture, ArbitraryOrder); using TestTypes = gtest_list_t< tuple_split_t<50, cross_product_t>>; INSTANTIATE_TYPED_TEST_CASE_P(Split1, TestGaussIntegrationFixture, TestTypes); using TestTypesTail = gtest_list_t< tuple_split_tail_t<50, cross_product_t>>; INSTANTIATE_TYPED_TEST_CASE_P(Split2, TestGaussIntegrationFixture, TestTypesTail); } diff --git a/test/test_fe_engine/test_fe_engine_precomputation.cc b/test/test_fe_engine/test_fe_engine_precomputation.cc index 46d178384..05668ed61 100644 --- a/test/test_fe_engine/test_fe_engine_precomputation.cc +++ b/test/test_fe_engine/test_fe_engine_precomputation.cc @@ -1,112 +1,111 @@ /** * @file test_fe_engine_precomputation.cc * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 - * @date last modification: Mon Jul 13 2015 + * @date last modification: Mon Feb 19 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "pybind11_akantu.hh" #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace py = pybind11; using namespace py::literals; template class TestFEMPyFixture : public TestFEMFixture { using parent = TestFEMFixture; public: void SetUp() override { parent::SetUp(); const auto & connectivities = this->mesh->getConnectivity(this->type); const auto & nodes = this->mesh->getNodes().begin(this->dim); coordinates = std::make_unique>( connectivities.size(), connectivities.getNbComponent() * this->dim); for (auto && tuple : zip(make_view(connectivities, connectivities.getNbComponent()), make_view(*coordinates, this->dim, connectivities.getNbComponent()))) { const auto & conn = std::get<0>(tuple); const auto & X = std::get<1>(tuple); for (auto s : arange(conn.size())) { Vector(X(s)) = Vector(nodes[conn(s)]); } } } void TearDown() override { parent::TearDown(); coordinates.reset(nullptr); } protected: std::unique_ptr> coordinates; }; TYPED_TEST_CASE(TestFEMPyFixture, types); TYPED_TEST(TestFEMPyFixture, Precompute) { SCOPED_TRACE(aka::to_string(this->type)); this->fem->initShapeFunctions(); const auto & N = this->fem->getShapeFunctions().getShapes(this->type); const auto & B = this->fem->getShapeFunctions().getShapesDerivatives(this->type); const auto & j = this->fem->getIntegrator().getJacobians(this->type); // Array ref_N(this->nb_quadrature_points_total, N.getNbComponent()); // Array ref_B(this->nb_quadrature_points_total, B.getNbComponent()); Array ref_j(this->nb_quadrature_points_total, j.getNbComponent()); auto ref_N(N); auto ref_B(B); py::module py_engine = py::module::import("py_engine"); auto py_shape = py_engine.attr("Shapes")(py::str(aka::to_string(this->type))); auto kwargs = py::dict( "N"_a = make_proxy(ref_N), "B"_a = make_proxy(ref_B), "j"_a = make_proxy(ref_j), "X"_a = make_proxy(*this->coordinates), "Q"_a = make_proxy(this->fem->getIntegrationPoints(this->type))); auto ret = py_shape.attr("precompute")(**kwargs); auto check = [&](auto & ref_A, auto & A, const auto & id) { SCOPED_TRACE(aka::to_string(this->type) + " " + id); for (auto && n : zip(make_view(ref_A, ref_A.getNbComponent()), make_view(A, A.getNbComponent()))) { auto diff = (std::get<0>(n) - std::get<1>(n)).template norm(); EXPECT_NEAR(0., diff, 1e-10); } }; check(ref_N, N, "N"); check(ref_B, B, "B"); check(ref_j, j, "j"); } diff --git a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc index 54804fda6..3df31e65e 100644 --- a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc +++ b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc @@ -1,153 +1,153 @@ /** - * @file test_fe_engine_precomputation.cc + * @file test_fe_engine_precomputation_bernoulli_2.cc * + * @author Lucas Frerot * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 - * @date last modification: Mon Jul 13 2015 + * @date last modification: Thu Jan 25 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; Matrix globalToLocalRotation(Real theta) { // clang-format off return {{ std::cos(theta), std::sin(theta), 0}, {-std::sin(theta), std::cos(theta), 0}, { 0, 0, 1}}; // clang-format on } Vector axialReference(Real xq) { return {(1. - xq) / 2, 0, 0, (1. + xq) / 2, 0, 0}; } Vector bendingReference(Real xq) { return {0, 1. / 4. * Math::pow<2>(xq - 1) * (xq + 2), 1. / 4. * Math::pow<2>(xq - 1) * (xq + 1), 0, 1. / 4. * Math::pow<2>(xq + 1) * (2 - xq), 1. / 4. * Math::pow<2>(xq + 1) * (xq - 1)}; } Vector bendingRotationReference(Real xq) { return {0, 3. / 4. * (xq * xq - 1), 1. / 4. * (3 * xq * xq - 2 * xq - 1), 0, 3. / 4. * (1 - xq * xq), 1. / 4. * (3 * xq * xq + 2 * xq - 1)}; } bool testBending(const Array & shape_functions, UInt shape_line_index, std::function(Real)> reference) { Real xq = -1. / std::sqrt(3.); // Testing values for bending rotations shapes on quadrature points for (auto && N : make_view(shape_functions, 3, 6)) { auto Nt = N.transpose(); Vector N_bending = Nt(shape_line_index); auto bending_reference = reference(xq); if (!Math::are_vector_equal(6, N_bending.storage(), bending_reference.storage())) return false; xq *= -1; } std::cout.flush(); return true; } int main(int argc, char * argv[]) { akantu::initialize(argc, argv); // debug::setDebugLevel(dblTest); constexpr ElementType type = _bernoulli_beam_2; UInt dim = ElementClass::getSpatialDimension(); Mesh mesh(dim); // creating nodes Vector node = {0, 0}; mesh.getNodes().push_back(node); node = {3. / 5., 4. / 5.}; mesh.getNodes().push_back(node); node = {2 * 3. / 5., 0}; mesh.getNodes().push_back(node); mesh.addConnectivityType(type); auto & connectivity = mesh.getConnectivity(type); // creating elements Vector elem = {0, 1}; connectivity.push_back(elem); elem = {1, 2}; connectivity.push_back(elem); elem = {0, 2}; connectivity.push_back(elem); using FE = FEEngineTemplate; using ShapeStruct = ShapeStructural<_ek_structural>; auto fem = std::make_unique(mesh, dim, "test_fem"); fem->initShapeFunctions(); auto & shape = dynamic_cast(fem->getShapeFunctions()); Array angles; angles.push_back(std::atan(4. / 3.)); angles.push_back(-std::atan(4. / 3.)); angles.push_back(0); /// Testing the rotation matrices for (auto && tuple : zip(make_view(shape.getRotations(type), 3, 3), angles)) { auto && rotation = std::get<0>(tuple); auto theta = std::get<1>(tuple); auto reference = globalToLocalRotation(theta); if (!Math::are_vector_equal(9, reference.storage(), rotation.storage())) return 1; } auto & shape_functions = shape.getShapes(type); if (!testBending(shape_functions, 0, axialReference)) return 1; // if (!testBending(shape_functions, 1, bendingReference)) // return 1; // if (!testBending(shape_functions, 2, bendingRotationReference)) // return 1; std::cout.flush(); finalize(); return 0; } diff --git a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc index 8f7828579..7b618d76c 100644 --- a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc +++ b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_3.cc @@ -1,106 +1,105 @@ /** - * @file test_fe_engine_precomputation.cc + * @file test_fe_engine_precomputation_bernoulli_3.cc * + * @author Lucas Frerot * @author Nicolas Richart * - * @date creation: Mon Jun 14 2010 - * @date last modification: Mon Jul 13 2015 + * @date creation: Wed Jan 24 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; /** * Reference: p. 285, example 5.7 - A First Course in the Finite Elements Method * Logan, 6th Edition, 2016 * ISBN-13: 978-1-305-63734-4 */ Matrix rotationReference() { return {{3. / 13, 4. / 13, 12. / 13}, {-4. / 5, 3. / 5, 0}, {-36. / 65, -48. / 65, 5. / 13}}; } int main(int argc, char * argv[]) { akantu::initialize(argc, argv); // debug::setDebugLevel(dblTest); constexpr ElementType type = _bernoulli_beam_3; UInt dim = ElementClass::getSpatialDimension(); Mesh mesh(dim); // Pushing nodes Vector node = {0, 0, 0}; mesh.getNodes().push_back(node); node = {3, 4, 12}; mesh.getNodes().push_back(node); // Pushing connectivity mesh.addConnectivityType(type); auto & connectivity = mesh.getConnectivity(type); Vector elem = {0, 1}; connectivity.push_back(elem); // Pushing normals auto & normals = mesh.registerData("extra_normal").alloc(0, dim, type, _not_ghost); Vector normal = {-36. / 65, -48. / 65, 5. / 13}; normals.push_back(normal); normals.push_back(normal); using FE = FEEngineTemplate; using ShapeStruct = ShapeStructural<_ek_structural>; auto fem = std::make_unique(mesh, dim, "test_fem"); fem->initShapeFunctions(); auto & shape = dynamic_cast(fem->getShapeFunctions()); Matrix rot_ref = rotationReference(); Matrix solution(6, 6); solution.block(rot_ref, 0, 0); solution.block(rot_ref, 3, 3); for (auto && rot : make_view(shape.getRotations(type), 6, 6)) { if (!Math::are_vector_equal(6 * 6, solution.storage(), rot.storage())) return 1; } /// TODO check shape functions and shape derivatives finalize(); return 0; } diff --git a/test/test_fe_engine/test_fe_engine_precomputation_structural.cc b/test/test_fe_engine/test_fe_engine_precomputation_structural.cc index 674f312e2..a396527ed 100644 --- a/test/test_fe_engine/test_fe_engine_precomputation_structural.cc +++ b/test/test_fe_engine/test_fe_engine_precomputation_structural.cc @@ -1,126 +1,126 @@ /** * @file test_fe_engine_precomputation_structural.cc * - * @author Lucas Frérot + * @author Lucas Frerot * - * @date creation: Fri Feb 26 2018 + * @date creation: Fri Jan 26 2018 + * @date last modification: Mon Feb 19 2018 * * @brief test of the structural precomputations * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_structural.hh" #include "test_fe_engine_structural_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ // Need a special fixture for the extra normal class TestBernoulliB3 : public TestFEMStructuralFixture> { using parent = TestFEMStructuralFixture>; public: /// Load the mesh and provide extra normal direction void readMesh(std::string filename) override { parent::readMesh(filename); auto & normals = this->mesh->registerData("extra_normal") .alloc(0, dim, type, _not_ghost); Vector normal = {-36. / 65, -48. / 65, 5. / 13}; normals.push_back(normal); } }; /* -------------------------------------------------------------------------- */ /// Type alias using TestBernoulliB2 = TestFEMStructuralFixture>; using TestDKT18 = TestFEMStructuralFixture>; /* -------------------------------------------------------------------------- */ /// Solution for 2D rotation matrices Matrix globalToLocalRotation(Real theta) { auto c = std::cos(theta); auto s = std::sin(theta); return {{c, s, 0}, {-s, c, 0}, {0, 0, 1}}; } /* -------------------------------------------------------------------------- */ TEST_F(TestBernoulliB2, PrecomputeRotations) { this->fem->initShapeFunctions(); using ShapeStruct = ShapeStructural<_ek_structural>; auto & shape = dynamic_cast(fem->getShapeFunctions()); auto & rot = shape.getRotations(type); Real a = std::atan(4. / 3); std::vector angles = {a, -a, 0}; Math::setTolerance(1e-15); for (auto && tuple : zip(make_view(rot, ndof, ndof), angles)) { auto rotation = std::get<0>(tuple); auto angle = std::get<1>(tuple); auto rotation_error = (rotation - globalToLocalRotation(angle)).norm(); EXPECT_NEAR(rotation_error, 0., Math::getTolerance()); } } /* -------------------------------------------------------------------------- */ TEST_F(TestBernoulliB3, PrecomputeRotations) { this->fem->initShapeFunctions(); using ShapeStruct = ShapeStructural<_ek_structural>; auto & shape = dynamic_cast(fem->getShapeFunctions()); auto & rot = shape.getRotations(type); Matrix ref = {{3. / 13, 4. / 13, 12. / 13}, {-4. / 5, 3. / 5, 0}, {-36. / 65, -48. / 65, 5. / 13}}; Matrix solution{ndof, ndof}; solution.block(ref, 0, 0); solution.block(ref, dim, dim); // The default tolerance is too much, really Math::setTolerance(1e-15); for (auto & rotation : make_view(rot, ndof, ndof)) { auto rotation_error = (rotation - solution).norm(); EXPECT_NEAR(rotation_error, 0., Math::getTolerance()); } } /* -------------------------------------------------------------------------- */ TEST_F(TestDKT18, PrecomputeRotations) { this->fem->initShapeFunctions(); using ShapeStruct = ShapeStructural<_ek_structural>; auto & shape = dynamic_cast(fem->getShapeFunctions()); auto & rot = shape.getRotations(type); for (auto & rotation : make_view(rot, ndof, ndof)) { std::cout << rotation << "\n"; } std::cout.flush(); } diff --git a/test/test_fe_engine/test_fe_engine_structural_fixture.hh b/test/test_fe_engine/test_fe_engine_structural_fixture.hh index 65c62c874..850d738f4 100644 --- a/test/test_fe_engine/test_fe_engine_structural_fixture.hh +++ b/test/test_fe_engine/test_fe_engine_structural_fixture.hh @@ -1,64 +1,64 @@ /** * @file test_fe_engine_structural_fixture.hh * - * @author Lucas Frérot + * @author Lucas Frerot * - * @date creation: Fri Feb 26 2018 + * @date creation: Fri Aug 20 2010 + * @date last modification: Fri Jan 26 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_io_msh_struct.hh" #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_FE_ENGINE_STRUCTURAL_FIXTURE_HH__ #define __AKANTU_TEST_FE_ENGINE_STRUCTURAL_FIXTURE_HH__ using namespace akantu; /// Base class for structural FEEngine tests with structural elements template class TestFEMStructuralFixture : public TestFEMBaseFixture { using parent = TestFEMBaseFixture; public: static const UInt ndof = ElementClass::getNbDegreeOfFreedom(); /// Need to tell the mesh to load structural elements void readMesh(std::string file_name) override { this->mesh->read(file_name, _miot_gmsh_struct); } }; template const UInt TestFEMStructuralFixture::ndof; // using types = gtest_list_t; // TYPED_TEST_CASE(TestFEMFixture, types); #endif /* __AKANTU_TEST_FE_ENGINE_STRUCTURAL_FIXTURE_HH__ */ diff --git a/test/test_fe_engine/test_gradient.cc b/test/test_fe_engine/test_gradient.cc index 521d548e1..9e8c13f27 100644 --- a/test/test_fe_engine/test_gradient.cc +++ b/test/test_fe_engine/test_gradient.cc @@ -1,106 +1,106 @@ /** * @file test_gradient.cc * * @author Nicolas Richart * @author Peter Spijker * * @date creation: Fri Sep 03 2010 - * @date last modification: Thu Oct 15 2015 + * @date last modification: Mon Feb 19 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * * @section DESCRIPTION * * This code is computing the gradient of a linear field and check that it gives * a constant result. It also compute the gradient the coordinates of the mesh * and check that it gives the identity * */ + /* -------------------------------------------------------------------------- */ #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { TYPED_TEST(TestFEMFixture, GradientPoly) { this->fem->initShapeFunctions(); Real alpha[2][3] = {{13, 23, 31}, {11, 7, 5}}; const auto dim = this->dim; const auto type = this->type; const auto & position = this->fem->getMesh().getNodes(); Array const_val(this->fem->getMesh().getNbNodes(), 2, "const_val"); for (auto && pair : zip(make_view(position, dim), make_view(const_val, 2))) { auto & pos = std::get<0>(pair); auto & const_ = std::get<1>(pair); const_.set(0.); for (UInt d = 0; d < dim; ++d) { const_(0) += alpha[0][d] * pos(d); const_(1) += alpha[1][d] * pos(d); } } /// compute the gradient Array grad_on_quad(this->nb_quadrature_points_total, 2 * dim, "grad_on_quad"); this->fem->gradientOnIntegrationPoints(const_val, grad_on_quad, 2, type); /// check the results for (auto && grad : make_view(grad_on_quad, 2, dim)) { for (UInt d = 0; d < dim; ++d) { EXPECT_NEAR(grad(0, d), alpha[0][d], 5e-13); EXPECT_NEAR(grad(1, d), alpha[1][d], 5e-13); } } } TYPED_TEST(TestFEMFixture, GradientPositions) { this->fem->initShapeFunctions(); const auto dim = this->dim; const auto type = this->type; UInt nb_quadrature_points = this->fem->getNbIntegrationPoints(type) * this->nb_element; Array grad_coord_on_quad(nb_quadrature_points, dim * dim, "grad_coord_on_quad"); const auto & position = this->mesh->getNodes(); this->fem->gradientOnIntegrationPoints(position, grad_coord_on_quad, dim, type); auto I = Matrix::eye(UInt(dim)); for (auto && grad : make_view(grad_coord_on_quad, dim, dim)) { auto diff = (I - grad).template norm(); EXPECT_NEAR(0., diff, 2e-14); } } } diff --git a/test/test_fe_engine/test_integrate.cc b/test/test_fe_engine/test_integrate.cc index 74ad310ad..a7b05be92 100644 --- a/test/test_fe_engine/test_integrate.cc +++ b/test/test_fe_engine/test_integrate.cc @@ -1,77 +1,77 @@ /** * @file test_integrate.cc * * @author Guillaume Anciaux * @author Nicolas Richart * @author Peter Spijker * * @date creation: Fri Sep 03 2010 - * @date last modification: Thu Oct 15 2015 + * @date last modification: Mon Feb 19 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ + /* -------------------------------------------------------------------------- */ #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { TYPED_TEST(TestFEMFixture, IntegrateConstant) { this->fem->initShapeFunctions(); const auto type = this->type; const auto & position = this->fem->getMesh().getNodes(); Array const_val(position.size(), 2, "const_val"); Array val_on_quad(this->nb_quadrature_points_total, 2, "val_on_quad"); Vector value{1, 2}; for (auto && const_ : make_view(const_val, 2)) { const_ = value; } // interpolate function on quadrature points this->fem->interpolateOnIntegrationPoints(const_val, val_on_quad, 2, type); // integrate function on elements Array int_val_on_elem(this->nb_element, 2, "int_val_on_elem"); this->fem->integrate(val_on_quad, int_val_on_elem, 2, type); // get global integration value Vector sum{0., 0.}; for (auto && int_ : make_view(int_val_on_elem, 2)) { sum += int_; } auto diff = (value - sum).template norm(); EXPECT_NEAR(0, diff, 1e-14); } } // namespace diff --git a/test/test_fe_engine/test_interpolate.cc b/test/test_fe_engine/test_interpolate.cc index 77e2b2d1a..30479fa03 100644 --- a/test/test_fe_engine/test_interpolate.cc +++ b/test/test_fe_engine/test_interpolate.cc @@ -1,73 +1,72 @@ /** * @file test_interpolate.cc * * @author Nicolas Richart * * @date creation: Fri Sep 03 2010 - * @date last modification: Thu Oct 15 2015 + * @date last modification: Tue Nov 14 2017 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { TYPED_TEST(TestFEMFixture, InterpolateConstant) { const auto type = this->type; const auto & position = this->fem->getMesh().getNodes(); Array const_val(position.size(), 2, "const_val"); Array val_on_quad(this->nb_quadrature_points_total, 2, "val_on_quad"); Vector value{1, 2}; for (auto && const_ : make_view(const_val, 2)) { const_ = value; } // interpolate function on quadrature points this->fem->interpolateOnIntegrationPoints(const_val, val_on_quad, 2, type); for (auto && int_ : make_view(val_on_quad, 2)) { auto diff = (value - int_).template norm(); EXPECT_NEAR(0, diff, 1e-14); } } // TYPED_TEST(TestFEMFixture, InterpolatePosition) { // const auto dim = this->dim; // const auto type = this->type; // const auto & position = this->fem->getMesh().getNodes(); // Array coord_on_quad(this->nb_quadrature_points_total, dim, // "coord_on_quad"); // this->fem->interpolateOnIntegrationPoints(position, coord_on_quad, dim, // type); // } } // namespace diff --git a/test/test_fe_engine/test_interpolate_bernoulli_beam_2.cc b/test/test_fe_engine/test_interpolate_bernoulli_beam_2.cc index 8772d060f..f81b4fde6 100644 --- a/test/test_fe_engine/test_interpolate_bernoulli_beam_2.cc +++ b/test/test_fe_engine/test_interpolate_bernoulli_beam_2.cc @@ -1,119 +1,118 @@ /** * @file test_interpolate_bernoulli_beam_2.cc * * @author Fabian Barras * * @date creation: Fri Jul 15 2011 - * @date last modification: Mon Jul 13 2015 + * @date last modification: Sat Jan 23 2016 * * @brief Test of the interpolation on the type _bernoulli_beam_2 * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "fe_engine.hh" #include "fe_engine_template.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "shape_linked.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main() { Mesh beams(2); /* -------------------------------------------------------------------------- */ // Defining the mesh Array & nodes = const_cast &>(beams.getNodes()); nodes.resize(4); beams.addConnectivityType(_bernoulli_beam_2); Array & connectivity = const_cast &>(beams.getConnectivity(_bernoulli_beam_2)); connectivity.resize(3); for (UInt i = 0; i < 4; ++i) { nodes(i, 0) = (i + 1) * 2; nodes(i, 1) = 1; } for (UInt i = 0; i < 3; ++i) { connectivity(i, 0) = i; connectivity(i, 1) = i + 1; } akantu::MeshIOMSH mesh_io; mesh_io.write("b_beam_2.msh", beams); /* -------------------------------------------------------------------------- */ // Interpolation FEEngineTemplate * fem = new FEEngineTemplate(beams, 2); fem->initShapeFunctions(); Array displ_on_nodes(4, 3); Array displ_on_quad(0, 3); for (UInt i = 0; i < 4; ++i) { displ_on_nodes(i, 0) = (i + 1) * 2; // Definition of the displacement displ_on_nodes(i, 1) = 0; displ_on_nodes(i, 2) = 0; } fem->getShapeFunctions().interpolateOnControlPoints<_bernoulli_beam_2>( displ_on_nodes, displ_on_quad, 3, _not_ghost, NULL, false, 0, 0, 0); fem->getShapeFunctions().interpolateOnControlPoints<_bernoulli_beam_2>( displ_on_nodes, displ_on_quad, 3, _not_ghost, NULL, false, 1, 1, 1); fem->getShapeFunctions().interpolateOnControlPoints<_bernoulli_beam_2>( displ_on_nodes, displ_on_quad, 3, _not_ghost, NULL, true, 2, 2, 1); fem->getShapeFunctions().interpolateOnControlPoints<_bernoulli_beam_2>( displ_on_nodes, displ_on_quad, 3, _not_ghost, NULL, false, 3, 2, 3); fem->getShapeFunctions().interpolateOnControlPoints<_bernoulli_beam_2>( displ_on_nodes, displ_on_quad, 3, _not_ghost, NULL, true, 4, 3, 3); Real * don = displ_on_nodes.storage(); Real * doq = displ_on_quad.storage(); std::ofstream my_file("out.txt"); my_file << don << std::endl; my_file << doq << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_fe_engine/test_inverse_map.cc b/test/test_fe_engine/test_inverse_map.cc index 32d62e4c9..cb3297c0b 100644 --- a/test/test_fe_engine/test_inverse_map.cc +++ b/test/test_fe_engine/test_inverse_map.cc @@ -1,74 +1,74 @@ /** * @file test_inverse_map.cc * * @author Guillaume Anciaux * * @date creation: Fri Sep 03 2010 - * @date last modification: Thu Oct 15 2015 + * @date last modification: Mon Feb 19 2018 * * @brief test of the fem class * * @section LICENSE * - * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) + * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ + /* -------------------------------------------------------------------------- */ #include "test_fe_engine_fixture.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; namespace { TYPED_TEST(TestFEMFixture, InverseMap) { this->fem->initShapeFunctions(); Matrix quad = GaussIntegrationElement::getQuadraturePoints(); const auto & position = this->fem->getMesh().getNodes(); /// get the quadrature points coordinates Array coord_on_quad(quad.cols() * this->nb_element, this->dim, "coord_on_quad"); this->fem->interpolateOnIntegrationPoints(position, coord_on_quad, this->dim, this->type); Vector natural_coords(this->dim); auto length = (this->upper - this->lower).template norm(); for (auto && enum_ : enumerate(make_view(coord_on_quad, this->dim, quad.cols()))) { auto el = std::get<0>(enum_); const auto & quads_coords = std::get<1>(enum_); for (auto q : arange(quad.cols())) { Vector quad_coord = quads_coords(q); Vector ref_quad_coord = quad(q); this->fem->inverseMap(quad_coord, el, this->type, natural_coords); auto dis_normalized = ref_quad_coord.distance(natural_coords) / length; EXPECT_NEAR(0., dis_normalized, 5e-12); } } } } // namespace diff --git a/test/test_fe_engine/test_mesh_boundary.cc b/test/test_fe_engine/test_mesh_boundary.cc index 9e34b4e99..496b78b84 100644 --- a/test/test_fe_engine/test_mesh_boundary.cc +++ b/test/test_fe_engine/test_mesh_boundary.cc @@ -1,64 +1,65 @@ /** * @file test_mesh_boundary.cc * * @author Dana Christen * * @date creation: Fri May 03 2013 - * @date last modification: Mon Jul 13 2015 + * @date last modification: Sun Aug 13 2017 * * @brief Thest the element groups * * @section LICENSE * - * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ + /* -------------------------------------------------------------------------- */ #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { UInt spatialDimension(3); akantu::initialize(argc, argv); Mesh mesh(spatialDimension, "mesh_names"); std::cout << "Loading the mesh." << std::endl; mesh.read("./cube_physical_names.msh"); std::cout << "Examining mesh:" << std::endl; // Inspection of the number of boundaries UInt nb_boundaries = mesh.getNbElementGroups(spatialDimension - 1); AKANTU_DEBUG_INFO(nb_boundaries << " boundaries advertised by Mesh."); if (nb_boundaries == 0) { std::cout << "No boundary detected!" << std::endl; return 1; } std::cout << (*dynamic_cast(&mesh)) << std::endl; akantu::finalize(); return 0; } diff --git a/test/test_fe_engine/test_mesh_data.cc b/test/test_fe_engine/test_mesh_data.cc index 561c98aec..959bb6da1 100644 --- a/test/test_fe_engine/test_mesh_data.cc +++ b/test/test_fe_engine/test_mesh_data.cc @@ -1,86 +1,86 @@ /** * @file test_mesh_data.cc * * @author Dana Christen * * @date creation: Fri May 03 2013 - * @date last modification: Mon Jul 13 2015 + * @date last modification: Wed Feb 03 2016 * * @brief Test of the MeshData class * * @section LICENSE * - * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ #include #include #define QUOTES(x) #x #define ADD_QUOTES(x) QUOTES(x) #define CAT(x, y) x##_##y #define CONCAT(x, y) CAT(x, y) //#define TYPE std::string //#define VALUE1 "abc" //#define VALUE2 "qwe" #define ELEMENT _triangle_6 #define NAME CONCAT(TYPE, data) /* -------------------------------------------------------------------------- */ using namespace akantu; using namespace std; int main() { std::cout << "Testing with type " << ADD_QUOTES(TYPE) << " and values " << ADD_QUOTES(VALUE1) << "," << ADD_QUOTES(VALUE2) << "..." << std::endl; MeshData mesh_data; ElementType elem_type = ELEMENT; const std::string name = ADD_QUOTES(NAME); Array & vec = mesh_data.getElementalDataArrayAlloc(name, elem_type); // XXX TO DELETE // vec.copy(mesh_data.getElementalDataArrayAlloc(name, elem_type)); TYPE value[2] = {VALUE1, VALUE2}; vec.push_back(value[0]); vec.push_back(value[1]); for (UInt i(0); i < 2; i++) { AKANTU_DEBUG_ASSERT(vec(i) == value[i], "The Array accessed through the " "getElementDataArray method does " "not contain the right value."); } std::cout << vec << std::endl; std::cout << mesh_data.getTypeCode(name) << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_geometry/CMakeLists.txt b/test/test_geometry/CMakeLists.txt index 8d8159008..32a503a45 100644 --- a/test/test_geometry/CMakeLists.txt +++ b/test/test_geometry/CMakeLists.txt @@ -1,55 +1,46 @@ #=============================================================================== # @file CMakeLists.txt # # @author Lucas Frerot # -# @date creation: Fri Sep 03 2010 -# @date last modification: Tue Jan 19 2016 +# @date creation: Fri Jun 19 2015 +# @date last modification: Sat Jan 23 2016 # # @brief configuration for solver tests # # @section LICENSE # -# Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de -# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des -# Solides) +# Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # -# Akantu 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. +# Akantu 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. # -# Akantu 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. +# Akantu 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 Akantu. If not, see . +# You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== register_test(test_geometry_predicates SOURCES test_geometry_predicates.cc PACKAGE CGAL ) register_test(test_geometry_intersection SOURCES test_geometry_intersection.cc FILES_TO_COPY test_geometry_triangle.msh PACKAGE CGAL ) register_test(test_segment_intersection_triangle_3 SOURCES test_segment_intersection_triangle_3.cc FILES_TO_COPY test_geometry_triangle.msh PACKAGE CGAL ) register_test(test_segment_intersection_tetrahedron_4 SOURCES test_segment_intersection_tetrahedron_4.cc FILES_TO_COPY test_geometry_tetrahedron.msh PACKAGE CGAL ) diff --git a/test/test_geometry/test_geometry_intersection.cc b/test/test_geometry/test_geometry_intersection.cc index 96b76a626..dc22d6500 100644 --- a/test/test_geometry/test_geometry_intersection.cc +++ b/test/test_geometry/test_geometry_intersection.cc @@ -1,131 +1,131 @@ /** * @file test_geometry_intersection.cc * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Feb 27 2015 - * @date last modification: Thu Jan 14 2016 + * @date last modification: Wed Jan 31 2018 * * @brief Tests the intersection module * * @section LICENSE * - * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory - * (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "geom_helper_functions.hh" #include "mesh_geom_factory.hh" #include "tree_type_helper.hh" #include "mesh_geom_common.hh" #include #include #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; typedef cgal::Cartesian K; typedef IntersectionTypeHelper, K>, K::Segment_3>::intersection_type result_type; typedef cgal::Spherical SK; typedef boost::variant> sk_inter_res; /*typedef CGAL::cpp11::result_of >)>::type sk_res;*/ typedef std::pair pair_type; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblWarning); Mesh mesh(2); mesh.read("test_geometry_triangle.msh"); MeshGeomFactory<2, _triangle_3, Triangle, K> factory(mesh); factory.constructData(); const TreeTypeHelper, K>::tree & tree = factory.getTree(); K::Point_3 a(0., 0.25, 0.), b(1., 0.25, 0.); K::Segment_3 line(a, b); K::Point_3 begin(a), intermediate(0.25, 0.25, 0.), end(0.75, 0.25, 0.); K::Segment_3 result_0(begin, intermediate), result_1(intermediate, end); std::list list_of_intersections; tree.all_intersections(line, std::back_inserter(list_of_intersections)); const result_type & intersection_0 = list_of_intersections.back(); const result_type & intersection_1 = list_of_intersections.front(); if (!intersection_0 || !intersection_1) return EXIT_FAILURE; /// *-> first is the intersection ; *->second is the primitive id if (const K::Segment_3 * segment = boost::get(&(intersection_0->first))) { if (!compareSegments(*segment, result_0)) { return EXIT_FAILURE; } } else return EXIT_FAILURE; if (const K::Segment_3 * segment = boost::get(&(intersection_1->first))) { if (!compareSegments(*segment, result_1)) { return EXIT_FAILURE; } } else return EXIT_FAILURE; SK::Sphere_3 sphere(SK::Point_3(0, 0, 0), 3.); SK::Segment_3 seg(SK::Point_3(0, 0, 0), SK::Point_3(2., 2., 2.)); SK::Line_arc_3 arc(seg); std::list s_results; CGAL::intersection(arc, sphere, std::back_inserter(s_results)); if (pair_type * pair = boost::get(&s_results.front())) { std::cout << "xi = " << to_double(pair->first.x()) << ", yi = " << to_double(pair->first.y()) << std::endl; if (!comparePoints(pair->first, SK::Circular_arc_point_3(1.0, 1.0, 1.0))) return EXIT_FAILURE; } else return EXIT_FAILURE; MeshGeomFactory<2, _triangle_3, Line_arc, SK> Sfactory(mesh); Sfactory.constructData(); finalize(); return EXIT_SUCCESS; } diff --git a/test/test_geometry/test_geometry_predicates.cc b/test/test_geometry/test_geometry_predicates.cc index 0a79bcb22..976ae4725 100644 --- a/test/test_geometry/test_geometry_predicates.cc +++ b/test/test_geometry/test_geometry_predicates.cc @@ -1,88 +1,88 @@ /** * @file test_geometry_predicates.cc * * @author Lucas Frerot * * @date creation: Fri Jan 04 2013 - * @date last modification: Thu Jan 14 2016 + * @date last modification: Wed Jan 31 2018 * * @brief Tests the geometry predicates * * @section LICENSE * - * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "geom_helper_functions.hh" #include "mesh_geom_common.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; typedef cgal::Cartesian K; typedef K::Point_3 Point; typedef K::Segment_3 Segment; int main(int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblWarning); Point a(0, 1, 0); Point b(0, 1, 1); Segment seg1(a, b); Segment seg2(b, a); if (!compareSegments(seg1, seg2)) return EXIT_FAILURE; // Testing sort + unique on list of segments std::vector> pair_list; pair_list.push_back(std::make_pair(seg1, 1)); pair_list.push_back(std::make_pair(seg2, 2)); segmentPairsLess sorter; std::sort(pair_list.begin(), pair_list.end(), sorter); std::vector>::iterator it = std::unique(pair_list.begin(), pair_list.end(), compareSegmentPairs); if (it - pair_list.begin() != 1) { std::cout << pair_list.size() << std::endl; return EXIT_FAILURE; } // Testing insertion in set std::set, segmentPairsLess> pair_set; pair_set.insert(pair_set.begin(), std::make_pair(seg1, 1)); pair_set.insert(pair_set.begin(), std::make_pair(seg2, 2)); if (pair_set.size() != 1) { std::cout << pair_set.size() << std::endl; return EXIT_FAILURE; } finalize(); return EXIT_SUCCESS; } diff --git a/test/test_geometry/test_segment_intersection_tetrahedron_4.cc b/test/test_geometry/test_segment_intersection_tetrahedron_4.cc index 36a3d54b2..1bee35e55 100644 --- a/test/test_geometry/test_segment_intersection_tetrahedron_4.cc +++ b/test/test_geometry/test_segment_intersection_tetrahedron_4.cc @@ -1,145 +1,145 @@ /** * @file test_segment_intersection_tetrahedron_4.cc * * @author Lucas Frerot * * @date creation: Fri Feb 27 2015 - * @date last modification: Thu Jan 14 2016 + * @date last modification: Wed Jan 31 2018 * * @brief Tests the intersection module with _tetrahedron_4 elements * * @section LICENSE * - * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory - * (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_segment_intersector.hh" #include "mesh_geom_common.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; typedef cgal::Cartesian K; typedef K::Point_3 Point; typedef K::Segment_3 Segment; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblError); Mesh mesh(3), interface_mesh(3, "interface_mesh"); mesh.read("test_geometry_tetrahedron.msh"); MeshSegmentIntersector<3, _tetrahedron_4> intersector(mesh, interface_mesh); intersector.constructData(); // Testing a segment going through the cube Point point(1., 1., 1.); Segment segment(CGAL::ORIGIN, point); intersector.computeIntersectionQuery(segment); std::cout << "number of seg_2 : " << interface_mesh.getNbElement(_segment_2) << std::endl; if (interface_mesh.getNbElement(_segment_2) != 2) return EXIT_FAILURE; Vector bary(2), bary1(2), bary2(2); Element test; test.type = _segment_2; test.element = 0; interface_mesh.getBarycenter(test, bary1); test.element = 1; interface_mesh.getBarycenter(test, bary2); Real first_bary[] = {1. / 6., 1. / 6., 1. / 6.}; Real second_bary[] = {2. / 3., 2. / 3., 2. / 3.}; // We don't know the order of the elements, so here we test permutations if (!((Math::are_vector_equal(3, bary1.storage(), first_bary) && Math::are_vector_equal(3, bary2.storage(), second_bary)) || (Math::are_vector_equal(3, bary1.storage(), second_bary) && Math::are_vector_equal(3, bary2.storage(), first_bary)))) return EXIT_FAILURE; // Testing a segment completely inside one element Point a(0.05, 0.05, 0.05), b(0.06, 0.06, 0.06); Segment inside_segment(a, b); intersector.computeIntersectionQuery(inside_segment); test.element = interface_mesh.getNbElement(_segment_2) - 1; interface_mesh.getBarycenter(test, bary); Real third_bary[] = {0.055, 0.055, 0.055}; if (!Math::are_vector_equal(3, bary.storage(), third_bary)) return EXIT_FAILURE; // Testing a segment whose end points are inside elements Point c(0.1, 0.1, 0.1), d(0.9, 0.9, 0.9); Segment crossing_segment(c, d); intersector.computeIntersectionQuery(crossing_segment); UInt el1 = interface_mesh.getNbElement(_segment_2) - 2; UInt el2 = el1 + 1; test.element = el1; interface_mesh.getBarycenter(test, bary1); test.element = el2; interface_mesh.getBarycenter(test, bary2); Real fourth_bary[] = {13. / 60., 13. / 60., 13. / 60.}; Real fifth_bary[] = {37. / 60., 37. / 60., 37. / 60.}; // We don't know the order of the elements, so here we test permutations if (!((Math::are_vector_equal(3, bary1.storage(), fourth_bary) && Math::are_vector_equal(3, bary2.storage(), fifth_bary)) || (Math::are_vector_equal(3, bary1.storage(), fifth_bary) && Math::are_vector_equal(3, bary2.storage(), fourth_bary)))) return EXIT_FAILURE; // Testing a segment along the edge of elements Point e(1, 0, 0), f(0, 1, 0); Segment edge_segment(e, f); UInt current_nb_elements = interface_mesh.getNbElement(_segment_2); intersector.computeIntersectionQuery(edge_segment); if (interface_mesh.getNbElement(_segment_2) != current_nb_elements + 1) return EXIT_FAILURE; test.element = interface_mesh.getNbElement(_segment_2) - 1; interface_mesh.getBarycenter(test, bary); Real sixth_bary[] = {0.5, 0.5, 0}; if (!Math::are_vector_equal(3, bary.storage(), sixth_bary)) return EXIT_FAILURE; return EXIT_SUCCESS; } diff --git a/test/test_geometry/test_segment_intersection_triangle_3.cc b/test/test_geometry/test_segment_intersection_triangle_3.cc index c837c417e..541c5ca65 100644 --- a/test/test_geometry/test_segment_intersection_triangle_3.cc +++ b/test/test_geometry/test_segment_intersection_triangle_3.cc @@ -1,138 +1,138 @@ /** * @file test_segment_intersection_triangle_3.cc * * @author Lucas Frerot * @author Clement Roux * * @date creation: Fri Feb 27 2015 - * @date last modification: Thu Jan 14 2016 + * @date last modification: Wed Jan 31 2018 * * @brief Tests the interface mesh generation * * @section LICENSE * - * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory - * (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 + * 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. * * Akantu 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 + * 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 Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "geom_helper_functions.hh" #include "mesh_geom_common.hh" #include "mesh_segment_intersector.hh" #include "mesh_sphere_intersector.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; typedef cgal::Cartesian K; typedef cgal::Spherical SK; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("", argc, argv); debug::setDebugLevel(dblError); Math::setTolerance(1e-10); Mesh mesh(2), interface_mesh(2, "interface_mesh"); mesh.read("test_geometry_triangle.msh"); MeshSegmentIntersector<2, _triangle_3> intersector(mesh, interface_mesh); intersector.constructData(); // Testing a segment going out of the mesh K::Point_3 a(0, 0.25, 0), b(1, 0.25, 0), c(0.25, 0, 0), d(0.25, 1, 0); K::Segment_3 h_interface(a, b), v_interface(c, d); std::list interface_list; interface_list.push_back(h_interface); interface_list.push_back(v_interface); intersector.computeIntersectionQueryList(interface_list); if (interface_mesh.getNbElement(_segment_2) != 4) return EXIT_FAILURE; Vector bary(2); Element test; test.element = 0; test.type = _segment_2; interface_mesh.getBarycenter(test, bary); Real first_bary[] = {0.125, 0.25}; if (!Math::are_vector_equal(2, bary.storage(), first_bary)) return EXIT_FAILURE; // Testing a segment completely inside an element K::Point_3 e(0.1, 0.33, 0), f(0.1, 0.67, 0); K::Segment_3 inside_segment(e, f); intersector.computeIntersectionQuery(inside_segment); test.element = interface_mesh.getNbElement(_segment_2) - 1; interface_mesh.getBarycenter(test, bary); Real second_bary[] = {0.1, 0.5}; if (!Math::are_vector_equal(2, bary.storage(), second_bary)) return EXIT_FAILURE; #if 0 // cgal::Spherical kernel testing the addition of nodes std::cout << "initial mesh size = " << mesh.getNodes().size() << " nodes" << std::endl; SK::Sphere_3 sphere(SK::Point_3(0, 1, 0), 0.2*0.2); SK::Sphere_3 sphere2(SK::Point_3(1, 0, 0), 0.4999999999); MeshSphereIntersector<2, _triangle_3> intersector_sphere(mesh); intersector_sphere.constructData(); std::list sphere_list; sphere_list.push_back(sphere); sphere_list.push_back(sphere2); intersector_sphere.computeIntersectionQueryList(sphere_list); std::cout << "final mesh size = " << mesh.getNodes().size() << std::endl; const Array new_node_triangle_3 = intersector_sphere.getNewNodePerElem(); const Array & nodes = mesh.getNodes(); std::cout << "New nodes :" << std::endl; std::cout << "node 5, x=" << nodes(4,0) << ", y=" << nodes(4,1) << std::endl; std::cout << "node 6, x=" << nodes(5,0) << ", y=" << nodes(5,1) << std::endl; std::cout << "node 7, x=" << nodes(6,0) << ", y=" << nodes(6,1) << std::endl; if ( (new_node_triangle_3(0,0) != 1) || (new_node_triangle_3(1,0) != 2)){ for(UInt k=0; k != new_node_triangle_3.size(); ++k){ std::cout << new_node_triangle_3(k,0) << " new nodes in element " << k << ", node(s): " << new_node_triangle_3(k,1) << ", " << new_node_triangle_3(k,3) << ", on segment(s):" << new_node_triangle_3(k,2) << ", " << new_node_triangle_3(k,4) << std::endl; } return EXIT_FAILURE; } #endif finalize(); return EXIT_SUCCESS; } diff --git a/test/test_gtest_main.cc b/test/test_gtest_main.cc index 32018c506..264c4d1a9 100644 --- a/test/test_gtest_main.cc +++ b/test/test_gtest_main.cc @@ -1,53 +1,83 @@ +/** + * @file test_gtest_main.cc + * + * @author Nicolas Richart + * + * @date creation: Thu Nov 09 2017 + * @date last modification: Tue Jan 16 2018 + * + * @brief Main function for gtest tests + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "communicator.hh" /* -------------------------------------------------------------------------- */ #include #if defined(AKANTU_TEST_USE_PYBIND11) #include namespace py = pybind11; #endif /* -------------------------------------------------------------------------- */ namespace { class AkaEnvironment : public ::testing::Environment { public: AkaEnvironment(int & argc, char **& argv) : argc(argc), argv(argv) {} // Override this to define how to set up the environment. void SetUp() override { ::akantu::initialize(argc, argv); #if defined(AKANTU_USE_PYBIND11) // py::initialize_interpreter(); #endif } // Override this to define how to tear down the environment. void TearDown() override { ::akantu::finalize(); #if defined(AKANTU_USE_PYBIND11) // py::finalize_interpreter(); #endif } protected: int & argc; char **& argv; }; } int main(int argc, char ** argv) { #if defined(AKANTU_TEST_USE_PYBIND11) py::scoped_interpreter guard{}; #endif ::testing::InitGoogleTest(&argc, argv); ::testing::AddGlobalTestEnvironment(new AkaEnvironment(argc, argv)); ::testing::TestEventListeners & listeners = ::testing::UnitTest::GetInstance()->listeners(); if (::akantu::Communicator::getStaticCommunicator().whoAmI() != 0) { delete listeners.Release(listeners.default_result_printer()); } return RUN_ALL_TESTS(); } diff --git a/test/test_gtest_utils.hh b/test/test_gtest_utils.hh index 995b8c8e3..526d2bcb2 100644 --- a/test/test_gtest_utils.hh +++ b/test/test_gtest_utils.hh @@ -1,203 +1,233 @@ +/** + * @file test_gtest_utils.hh + * + * @author Nicolas Richart + * + * @date creation: Tue Nov 14 2017 + * @date last modification: Wed Feb 21 2018 + * + * @brief Utils to help write tests + * + * @section LICENSE + * + * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu 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. + * + * Akantu 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 Akantu. If not, see . + * + */ + /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_iterators.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TEST_GTEST_UTILS_HH__ #define __AKANTU_TEST_GTEST_UTILS_HH__ namespace { /* -------------------------------------------------------------------------- */ template <::akantu::ElementType t> using element_type_t = std::integral_constant<::akantu::ElementType, t>; /* -------------------------------------------------------------------------- */ template struct gtest_list {}; template struct gtest_list> { using type = ::testing::Types; }; template using gtest_list_t = typename gtest_list::type; /* -------------------------------------------------------------------------- */ template struct tuple_concat {}; template struct tuple_concat, std::tuple> { using type = std::tuple; }; template using tuple_concat_t = typename tuple_concat::type; /* -------------------------------------------------------------------------- */ template