diff --git a/CMakeLists.txt b/CMakeLists.txt index 5df76930b..320e5d31c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,263 +1,264 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 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 . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # _ ,-, _ # ,--, /: :\/': :`\/: :\ # |`; ' `,' `.; `: | _ _ # | | | ' | |. | | | | # | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _ # | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | | # \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| | # `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_| # `==._ .. . /' # `-::-' #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.6) project(AKANTU) enable_language(CXX) #=============================================================================== # Misc. #=============================================================================== set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") INCLUDE(${AKANTU_CMAKE_DIR}/AkantuMacros.cmake) #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 0) set(AKANTU_MINOR_VERSION 1) set(AKANTU_BUILD_VERSION 0) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}" ) # Append the library version information to the library target properties if(NOT AKANTU_NO_LIBRARY_VERSION) set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES} VERSION "${AKANTU_VERSION}" SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}" ) endif(NOT AKANTU_NO_LIBRARY_VERSION) #=============================================================================== # Options #=============================================================================== option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) # compilation switch option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF) # features add_optional_package(IOHelper "Add IOHelper support in akantu" OFF) add_optional_package(MPI "Add MPI support in akantu" OFF) add_optional_package(Scotch "Add Scotch support in akantu" OFF) add_optional_package(Mumps "Add Mumps support in akantu" OFF) add_optional_package(BLAS "Use BLAS for arithmetic operations" OFF) add_optional_package(EPSN "Use EPSN streering environment" OFF) # Debug if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) #=========================================================================== # Dependencies #=========================================================================== find_package(Boost REQUIRED) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIRS}) endif() #check_for_isnan(AKANTU_ISNAN) #message("ISNAN = "${AKANTU_ISNAN}) #=============================================================================== # Library #=============================================================================== set(AKANTU_COMMON_SRC common/aka_common.cc common/aka_error.cc common/aka_extern.cc common/aka_static_memory.cc common/aka_memory.cc common/aka_vector.cc common/aka_math.cc fem/shape_lagrange.cc fem/integrator_gauss.cc fem/mesh.cc fem/fem.cc fem/element_class.cc fem/fem_template.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_boundary.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/material.cc model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_damage.cc mesh_utils/mesh_io.cc mesh_utils/mesh_pbc.cc mesh_utils/mesh_io/mesh_io_msh.cc +# mesh_utils/mesh_io/mesh_io_diana.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc solver/sparse_matrix.cc solver/solver.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) set(AKANTU_CONTACT_SRC model/solid_mechanics/contact.cc model/solid_mechanics/contact_search.cc model/solid_mechanics/contact_neighbor_structure.cc model/solid_mechanics/contact/contact_2d_explicit.cc model/solid_mechanics/contact/contact_search_2d_explicit.cc model/solid_mechanics/contact/regular_grid_neighbor_structure.cc model/solid_mechanics/contact/contact_search_explicit.cc model/solid_mechanics/contact/contact_3d_explicit.cc model/solid_mechanics/contact/grid_2d_neighbor_structure.cc model/solid_mechanics/contact/contact_rigid.cc model/solid_mechanics/contact/friction_coefficient.cc model/solid_mechanics/contact/friction_coefficient/unique_constant_fric_coef.cc model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef.cc model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef/ruina_slowness_fric_coef.cc ) if(AKANTU_BUILD_CONTACT) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} ${AKANTU_CONTACT_SRC}) endif() if(AKANTU_MPI_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} synchronizer/static_communicator_mpi.cc ) endif() if(AKANTU_SCOTCH_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} mesh_utils/mesh_partition/mesh_partition_scotch.cc ) endif() if(AKANTU_MUMPS_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} solver/solver_mumps.cc ) endif() set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integration_scheme model/solid_mechanics model/solid_mechanics/materials model/solid_mechanics/contact model/solid_mechanics/contact/friction_coefficient model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef synchronizer/ solver/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=========================================================================== # Documentation #=========================================================================== if(AKANTU_DOCUMENTATION) find_package(Doxygen REQUIRED) if(DOXYGEN_FOUND) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) add_subdirectory(doc/doxygen) endif(DOXYGEN_FOUND) endif(AKANTU_DOCUMENTATION) #=============================================================================== # Tests #=============================================================================== ENABLE_TESTING() INCLUDE(CTest) INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) #=============================================================================== # Examples #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif(AKANTU_EXAMPLES) diff --git a/common/aka_math.cc b/common/aka_math.cc index a278e37ce..788eccf8d 100644 --- a/common/aka_math.cc +++ b/common/aka_math.cc @@ -1,209 +1,211 @@ /** * @file aka_math.cc * @author Nicolas Richart * @date Wed Jul 28 12:13:46 2010 * * @brief Implementation of the math toolbox * * @section LICENSE * * Copyright (©) 2010-2011 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_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void Math::matrix_vector(UInt m, UInt n, const Vector & A, const Vector & x, - Vector & y) { + Vector & y, + Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.getSize() == x.getSize(), "The vector A(" << A.getID() << ") and the vector x(" << x.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * n, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT(x.getNbComponent() == n, "The vector x(" << x.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT(y.getNbComponent() == n, "The vector y(" << y.getID() << ") do not the good number of component."); UInt nb_element = A.getSize(); UInt offset_A = A.getNbComponent(); UInt offset_x = x.getNbComponent(); y.resize(nb_element); Real * A_val = A.values; Real * x_val = x.values; Real * y_val = y.values; for (UInt el = 0; el < nb_element; ++el) { - matrix_vector(m, n, A_val, x_val, y_val); + matrix_vector(m, n, A_val, x_val, y_val, alpha); A_val += offset_A; x_val += offset_x; y_val += offset_x; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::matrix_matrix(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, - Vector & C) { + Vector & C, + Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.getSize() == B.getSize(), "The vector A(" << A.getID() << ") and the vector B(" << B.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * k, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT(B.getNbComponent() == k * n , "The vector B(" << B.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT(C.getNbComponent() == m * n, "The vector C(" << C.getID() << ") do not the good number of component."); UInt nb_element = A.getSize(); UInt offset_A = A.getNbComponent(); UInt offset_B = B.getNbComponent(); UInt offset_C = C.getNbComponent(); C.resize(nb_element); Real * A_val = A.values; Real * B_val = B.values; Real * C_val = C.values; for (UInt el = 0; el < nb_element; ++el) { - matrix_matrix(m, n, k, A_val, B_val, C_val); + matrix_matrix(m, n, k, A_val, B_val, C_val, alpha); A_val += offset_A; B_val += offset_B; C_val += offset_C; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, - Vector & C) { + Vector & C, + Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.getSize() == B.getSize(), "The vector A(" << A.getID() << ") and the vector B(" << B.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * k, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT(B.getNbComponent() == k * n , "The vector B(" << B.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT(C.getNbComponent() == m * n, "The vector C(" << C.getID() << ") do not the good number of component."); UInt nb_element = A.getSize(); UInt offset_A = A.getNbComponent(); UInt offset_B = B.getNbComponent(); UInt offset_C = C.getNbComponent(); C.resize(nb_element); Real * A_val = A.values; Real * B_val = B.values; Real * C_val = C.values; for (UInt el = 0; el < nb_element; ++el) { - matrix_matrixt(m, n, k, A_val, B_val, C_val); + matrix_matrixt(m, n, k, A_val, B_val, C_val, alpha); A_val += offset_A; B_val += offset_B; C_val += offset_C; } AKANTU_DEBUG_OUT(); } -void Math::matrix33_eigenvalues(Real * A, Real *Adiag ){ +/* -------------------------------------------------------------------------- */ +void Math::matrix33_eigenvalues(Real * A, Real *Adiag) { ///d = determinant of Matrix A Real d = det3(A); ///b = trace of Matrix A Real b = A[0]+A[4]+A[8]; - + Real a = -1 ; /// c = 0.5*(trace(M^2)-trace(M)^2) Real c = A[3]*A[1] + A[2]*A[6] + A[5]*A[7] - A[0]*A[4] - A[0]*A[8] - A[4]*A[8]; /// Define x, y, z Real x = ((3*c/a) - ((b*b)/(a*a)))/3; Real y=((2*(b*b*b)/(a*a*a)) - (9*b*c/(a*a)) + (27*d/a))/27; Real z = (y*y)/4 + (x*x*x)/27; /// Define I, j, k, m, n, p (so equations are not so cluttered) Real i = sqrt(y*y/4 - z); Real j = -pow(i,1./3.); Real k; if (fabs(i)<1e-12) k = 0; else k = acos(-(y/(2*i))); - + Real m = cos(k/3); Real n = sqrt(3)*sin(k/3); Real p = b/(3*a); - Adiag[0]=-(2*j*m + p);; Adiag[1]=-(-j *(m + n) + p); Adiag[2]=-(-j * (m - n) + p); - } __END_AKANTU__ diff --git a/common/aka_math.hh b/common/aka_math.hh index 8f7515b5a..4c5130013 100644 --- a/common/aka_math.hh +++ b/common/aka_math.hh @@ -1,193 +1,210 @@ /** * @file aka_math.hh * @author Nicolas Richart * @date Wed Jul 28 11:51:56 2010 * * @brief mathematical operations * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_MATH_H__ #define __AKANTU_AKA_MATH_H__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_vector.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Math { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Matrix algebra */ /* ------------------------------------------------------------------------ */ /// @f$ y = A*x @f$ static void matrix_vector(UInt m, UInt n, const Vector & A, const Vector & x, - Vector & y); + Vector & y, Real alpha = 1.); /// @f$ y = A*x @f$ static inline void matrix_vector(UInt m, UInt n, const Real * A, const Real * x, - Real * y); + Real * y, Real alpha = 1.); + + /// @f$ y = A^t*x @f$ + static inline void matrixt_vector(UInt m, UInt n, + const Real * A, + const Real * x, + Real * y, Real alpha = 1.); /// @f$ C = A*B @f$ static void matrix_matrix(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, - Vector & C); + Vector & C, Real alpha = 1.); /// @f$ C = A*B^t @f$ static void matrix_matrixt(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, - Vector & C); + Vector & C, Real alpha = 1.); /// @f$ C = A*B @f$ static inline void matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C); + Real * C, Real alpha = 1.); /// @f$ C = A^t*B @f$ static inline void matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C); + Real * C, Real alpha = 1.); /// @f$ C = A*B^t @f$ static inline void matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C); + Real * C, Real alpha = 1.); /// @f$ C = A^t*B^t @f$ static inline void matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C); + Real * C, Real alpha = 1.); + + template + static inline void matMul(UInt m, UInt n, UInt k, + Real alpha, const Real * A, const Real * B, + Real beta, Real * C); + + template + static inline void matVectMul(UInt m, UInt n, + Real alpha, const Real * A, const Real * x, + Real beta, Real * y); + static void matrix33_eigenvalues(Real * A, Real * Adiag); /// determinent of a 3x3 matrix static inline Real det3(const Real * mat); /// determinent of a 2x2 matrix static inline Real det2(const Real * mat); /// inverse a 3x3 matrix static inline void inv3(const Real * mat, Real * inv); /// inverse a 2x2 matrix static inline void inv2(const Real * mat, Real * inv); /// vector cross product static inline void vectorProduct3(const Real * v1, const Real * v2, Real * res); /// compute normal a normal to a vector static inline void normal2(const Real * v1, Real * res); /// compute normal a normal to a vector static inline void normal3(const Real * v1,const Real * v2, Real * res); /// normalize a vector static inline void normalize2(Real * v); /// normalize a vector static inline void normalize3(Real * v); /// return norm of a 2-vector static inline Real norm2(const Real * v); /// return norm of a 3-vector static inline Real norm3(const Real * v); /// return the dot product between 2 vectors in 2d static inline Real vectorDot2(const Real * v1, const Real * v2); /// return the dot product between 2 vectors in 3d static inline Real vectorDot3(const Real * v1, const Real * v2); /* ------------------------------------------------------------------------ */ /* Geometry */ /* ------------------------------------------------------------------------ */ /// distance in 2D between x and y static inline Real distance_2d(const Real * x, const Real * y); /// distance in 3D between x and y static inline Real distance_3d(const Real * x, const Real * y); /// radius of the in-circle of a triangle static inline Real triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3); /// radius of the in-circle of a tetrahedron static inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4); /// volume of a tetrahedron static inline Real tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4); /// compute the barycenter of n points static inline void barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter); /// vector between x and y static inline void vector_2d(const Real * x, const Real * y, Real * vec); /// vector pointing from x to y in 3 spatial dimension static inline void vector_3d(const Real * x, const Real * y, Real * vec); /// test if two scalar are equal within a given tolerance static inline bool are_float_equal(Real x, Real y); /// tolerance for functions that need one static Real tolerance; /// test if a real is a NaN static inline bool isnan(Real x); }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "aka_math_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_AKA_MATH_H__ */ diff --git a/common/aka_math_inline_impl.cc b/common/aka_math_inline_impl.cc index 6524f8ed1..b4d99e485 100644 --- a/common/aka_math_inline_impl.cc +++ b/common/aka_math_inline_impl.cc @@ -1,396 +1,465 @@ /** * @file aka_math_inline_impl.cc * @author Nicolas Richart * @date Wed Jul 28 13:20:35 2010 * * @brief Implementation of the inline functions of the math toolkit * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ #ifdef AKANTU_USE_BLAS # ifndef AKANTU_USE_BLAS_MKL # include # else // AKANTU_USE_BLAS_MKL # include # endif //AKANTU_USE_BLAS_MKL #endif //AKANTU_USE_BLAS /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt m, UInt n, const Real * A, const Real * x, - Real * y) { + Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, - m, n, 1, A, n, x, 1, 0, y, 1); + m, n, alpha, A, n, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * n; for (UInt j = 0; j < n; ++j) { y[i] += A[A_i + j] * x[j]; } + y[i] *= alpha; + } +#endif +} + + +/* -------------------------------------------------------------------------- */ +inline void Math::matrixt_vector(UInt m, UInt n, + const Real * A, + const Real * x, + Real * y, Real alpha) { +#ifdef AKANTU_USE_BLAS + /// y = alpha*op(A)*x + beta*y + cblas_dgemv(CblasRowMajor, CblasNoTrans, + m, n, alpha, A, m, x, 1, 0, y, 1); +#else + memset(y, 0, m*sizeof(Real)); + for (UInt i = 0; i < m; ++i) { + Real xt = alpha * x[i]; + for (UInt j = 0; j < n; ++j) { + y[j] += A[i + j * m] * xt; + } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C) { + Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, - 1, + alpha, A, k, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[l * n + j]; } + C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C) { + Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, - 1, + alpha, A, m, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { // UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[l * n + j]; } + C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C) { + Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, - 1, + alpha, A, k, B, k, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[B_j + l]; } + C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, - Real * C) { + Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, m, n, k, - 1, + alpha, A, m, B, k, 0, C, n); #else memset(C, 0, m * n * sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[B_j + l]; } + C[C_i + j] *= alpha; } } #endif } +/* -------------------------------------------------------------------------- */ +template +inline void Math::matMul(UInt m, UInt n, UInt k, + Real alpha, const Real * A, const Real * B, + Real beta, Real * C) { + // #ifndef AKANTU_USE_BLAS + if(tr_A) { + if(tr_B) matrixt_matrixt(m, n, k, A, B, C, alpha); + else matrixt_matrix(m, n, k, A, B, C, alpha); + } else { + if(tr_B) matrix_matrixt(m, n, k, A, B, C, alpha); + else matrix_matrix(m, n, k, A, B, C, alpha); + } + // #else + // static CBLAS_TRANSPOSE transpose[2] = {CblasNoTrans, CblasTrans}; + // cblas_dgemm(CblasRowMajor, transpose[tr_A], transpose[tr_B], + // m, n, k, + // alpha, A, m, B, k, + // beta, C, n); + // #endif +} + +/* -------------------------------------------------------------------------- */ +template +inline void Math::matVectMul(UInt m, UInt n, + Real alpha, const Real * A, const Real * x, + Real beta, Real * y) { + //#ifndef AKANTU_USE_BLAS + if(tr_A) { + matrixt_vector(m, n, A, x, y, alpha); + } else { + matrix_vector(m, n, A, x, y, alpha); + } + // #else + // static CBLAS_TRANSPOSE transpose[2] = {CblasNoTrans, CblasTrans}; + // cblas_dgemv(CblasRowMajor, transpose[tr_A], + // m, n, + // alpha, A, m, x, 1, + // beta, y, 1); + // #endif +} + + /* -------------------------------------------------------------------------- */ inline Real Math::det2(const Real * mat) { return mat[0]*mat[3] - mat[1]*mat[2]; } /* -------------------------------------------------------------------------- */ inline Real Math::det3(const Real * mat) { return mat[0]*(mat[4]*mat[8]-mat[7]*mat[5]) - mat[3]*(mat[1]*mat[8]-mat[7]*mat[2]) + mat[6]*(mat[1]*mat[5]-mat[4]*mat[2]); } /* -------------------------------------------------------------------------- */ inline void Math::normal2(const Real * vec,Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; Math::normalize2(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normal3(const Real * vec1,const Real * vec2,Real * normal) { Math::vectorProduct3(vec1,vec2,normal); Math::normalize3(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normalize2(Real * vec) { Real norm = Math::norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* -------------------------------------------------------------------------- */ inline void Math::normalize3(Real * vec) { Real norm = Math::norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* -------------------------------------------------------------------------- */ inline Real Math::norm2(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm3(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); } /* -------------------------------------------------------------------------- */ inline void Math::inv2(const Real * mat,Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::inv3(const Real * mat,Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4]*mat[8] - mat[7]*mat[5])/det_mat; inv[1] = (mat[2]*mat[7] - mat[8]*mat[1])/det_mat; inv[2] = (mat[1]*mat[5] - mat[4]*mat[2])/det_mat; inv[3] = (mat[5]*mat[6] - mat[8]*mat[3])/det_mat; inv[4] = (mat[0]*mat[8] - mat[6]*mat[2])/det_mat; inv[5] = (mat[2]*mat[3] - mat[5]*mat[0])/det_mat; inv[6] = (mat[3]*mat[7] - mat[6]*mat[4])/det_mat; inv[7] = (mat[1]*mat[6] - mat[7]*mat[0])/det_mat; inv[8] = (mat[0]*mat[4] - mat[3]*mat[1])/det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot2(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot3(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_2d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1])); } /* -------------------------------------------------------------------------- */ inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\ * s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; a = distance_2d(coord1, coord2); b = distance_2d(coord2, coord3); c = distance_2d(coord1, coord3); Real s; s = (a + b + c) * 0.5; return sqrt((s - a) * (s - b) * (s - c) / s); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_3d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1]) + (y[2] - x[2])*(y[2] - x[2]) ); } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9], vol; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real l12, l13, l14, l23, l24, l34; l12 = distance_3d(coord1, coord2); l13 = distance_3d(coord1, coord3); l14 = distance_3d(coord1, coord4); l23 = distance_3d(coord2, coord3); l24 = distance_3d(coord2, coord4); l34 = distance_3d(coord3, coord4); Real s1, s2, s3, s4; s1 = (l12 + l23 + l13) * 0.5; s1 = sqrt(s1*(s1-l12)*(s1-l23)*(s1-l13)); s2 = (l12 + l24 + l14) * 0.5; s2 = sqrt(s2*(s2-l12)*(s2-l24)*(s2-l14)); s3 = (l23 + l34 + l24) * 0.5; s3 = sqrt(s3*(s3-l23)*(s3-l34)*(s3-l24)); s4 = (l13 + l34 + l24) * 0.5; s4 = sqrt(s4*(s4-l13)*(s4-l34)*(s4-l14)); Real volume = Math::tetrahedron_volume(coord1,coord2,coord3,coord4); return 3*volume/(s1+s2+s3+s4); } /* -------------------------------------------------------------------------- */ inline void Math::barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { memset(barycenter, 0, spatial_dimension * sizeof(Real)); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real) nb_points; } } } /* -------------------------------------------------------------------------- */ inline void Math::vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; } /* -------------------------------------------------------------------------- */ inline void Math::vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; res[2] = y[2]-x[2]; } /* -------------------------------------------------------------------------- */ inline bool Math::are_float_equal(const Real x, const Real y){ return (fabs( x - y) < tolerance); } /* -------------------------------------------------------------------------- */ inline bool Math::isnan(Real x) { return ::isnan(x); } diff --git a/common/aka_types.hh b/common/aka_types.hh index 1dcae032a..c4870cd89 100644 --- a/common/aka_types.hh +++ b/common/aka_types.hh @@ -1,311 +1,382 @@ /** * @file aka_types.hh * @author Nicolas Richart * @date Wed Feb 16 20:28:13 2011 * * @brief description of the "simple" types * * @section LICENSE * * Copyright (©) 2010-2011 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_vector.hh" #ifdef AKANTU_USE_BLAS # ifndef AKANTU_USE_CBLAS_MKL # include # else // AKANTU_USE_CBLAS_MKL # include # endif //AKANTU_USE_CBLAS_MKL #endif //AKANTU_USE_BLAS /* -------------------------------------------------------------------------- */ #ifndef __INTEL_COMPILER #include #endif /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_TYPES_HH__ #define __AKANTU_AKA_TYPES_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* maps */ /* -------------------------------------------------------------------------- */ #ifndef __INTEL_COMPILER template -struct unordered_map { - typedef typename std::tr1::unordered_map type; -}; +struct unordered_map { typedef typename std::tr1::unordered_map type; }; #else template -struct unordered_map { - typedef typename std::map type; -}; +struct unordered_map { typedef typename std::map type; }; #endif -/* -------------------------------------------------------------------------- */ -/* matrix */ -/* -------------------------------------------------------------------------- */ +namespace types { -template -class TMatrix { -public: - TMatrix() : values(new T[n*m]), wrapped(false) {}; + /* ------------------------------------------------------------------------ */ + /* Matrix */ + /* ------------------------------------------------------------------------ */ + class Matrix { + public: + Matrix() : m(0), n(0), values(NULL), wrapped(true) {}; - TMatrix(T def) : values(new T[n*m]), wrapped(false) { - std::fill_n(values, n*m, def); - }; + Matrix(UInt m, UInt n, Real def = 0) : m(m), n(n), values(new Real[n*m]), wrapped(false) { + std::fill_n(values, n*m, def); + }; - TMatrix(T *data) : values(data), wrapped(true) {}; + Matrix(Real* data, UInt m, UInt n) : m(m), n(n), values(data), wrapped(true) {}; - TMatrix(TMatrix & src) { - values = src.values; - wrapped = src.wrapped; - src.wrapped = true; - }; + Matrix(Matrix & src) { + values = src.values; + wrapped = src.wrapped; + src.wrapped = true; + }; - virtual ~TMatrix() { if(!wrapped) delete [] values; }; + virtual ~Matrix() { if(!wrapped) delete [] values; }; - static UInt size() { return n*m; }; + UInt size() const { return n*m; }; + UInt rows() const { return m; }; + UInt cols() const { return n; }; + Real * storage() const { return values; }; - inline T & operator()(UInt i, UInt j) { return *(values + i*m + j); }; - inline const T & operator()(UInt i, UInt j) const { return *(values + i*m + j); }; - inline T & operator[](UInt idx) { return *(values + idx); }; - inline const T & operator[](UInt idx) const { return *(values + idx); }; + inline Real& operator()(UInt i, UInt j) { return *(values + i*m + j); }; + inline const Real& operator()(UInt i, UInt j) const { return *(values + i*m + j); }; + inline Real& operator[](UInt idx) { return *(values + idx); }; + inline const Real& operator[](UInt idx) const { return *(values + idx); }; - inline TMatrix & operator=(const TMatrix & mat) { - if(this != &mat) { - if(values == NULL) { this->values = new T[n * m]; this->wrapped = false; } - memcpy(this->values, &mat[0], size()*sizeof(T)); - } - return *this; - }; + inline Matrix & operator=(const Matrix & mat) { + if(this != &mat) { + if(values != NULL) { + memcpy(this->values, mat.values, m*n*sizeof(Real)); + } else { + this->m = mat.m; + this->n = mat.n; - inline void clear() { memset(values, 0, m * n * sizeof(T)); }; + this->values = mat.values; + const_cast(mat).wrapped = true; + this->wrapped = false; + } + } + return *this; + }; + + inline Matrix operator* (const Matrix & B) { + // UInt m = this->m, n = B.n, k = this->n; + // Matrix C(m,n); + // C.clear(); + + // UInt A_i = 0, C_i = 0; + // for (UInt i = 0; i < m; ++i, A_i += k, C_i += n) + // for (UInt j = 0; j < n; ++j) + // for (UInt l = 0; l < k; ++l) + // C[C_i + j] += (*this)[A_i + l] * B[l * n + j]; + Matrix C(m,n); + C.mul(*this, B); + + return C; + }; + + inline void mul(const Matrix & A, const Matrix & B) { + UInt m = A.m, n = B.n, k = A.n; +#ifdef AKANTU_USE_BLAS + /// C := alpha*op(A)*op(B) + beta*C + cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, + m, n, k, + 1, A.values, k, B.values, n, + 0, this->values, n); +#else + this->clear(); + + UInt A_i = 0, C_i = 0; + for (UInt i = 0; i < m; ++i, A_i += k, C_i += n) + for (UInt j = 0; j < n; ++j) + for (UInt l = 0; l < k; ++l) + (*this)[C_i + j] += A[A_i + l] * B[l * n + j]; +#endif + }; - /// function to print the containt of the class - virtual void printself(std::ostream & stream, int indent = 0) const { - std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - stream << space << typeid(*this).name() << "," << n << "," << m <<"> :" << std::endl; - for (UInt i = 0; i < m; ++i) { - stream << space << space << "| "; - for (UInt j = 0; j < n; ++j) { - stream << values[i*n +j] << " "; + inline void clear() { memset(values, 0, m * n * sizeof(Real)); }; + + /// function to print the containt of the class + virtual void printself(std::ostream & stream, int indent = 0) const { + std::string space; + for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + + stream << space << typeid(*this).name() << "<" << n << "," << m <<"> :" << std::endl; + for (UInt i = 0; i < m; ++i) { + stream << space << indent << "| "; + for (UInt j = 0; j < n; ++j) { + stream << values[i*n +j] << " "; + } + stream << "|" << std::endl; } - stream << "|" << std::endl; - } - }; + }; - // T* &data() { return values; }; + // Real* &data() { return values; }; + // bool &is_wrapped() { return wrapped; }; + friend class ::akantu::Vector; - friend class Vector; + protected: + UInt m; + UInt n; + Real* values; + bool wrapped; + }; -protected: - T * values; - bool wrapped; -}; + /* ------------------------------------------------------------------------ */ + /* Vector */ + /* ------------------------------------------------------------------------ */ + class Vector { + public: + Vector() : n(0), values(NULL), wrapped(true) {}; -/* -------------------------------------------------------------------------- */ + Vector(UInt n, Real def = 0) : n(n), values(new Real[n]), wrapped(false) { + std::fill_n(values, n, def); + }; -template -class RealTMatrix : public TMatrix { -public: - RealTMatrix() : TMatrix() {}; - RealTMatrix(Real def) : TMatrix(def) {}; - RealTMatrix(Real * data) : TMatrix(data) {}; - RealTMatrix(RealTMatrix & mat) : TMatrix(mat) {}; -}; + Vector(Real* data, UInt n) : n(n), values(data), wrapped(true) {}; + Vector(Vector & src) { + values = src.values; + wrapped = src.wrapped; + src.wrapped = true; + }; -template -inline RealTMatrix operator* (const RealTMatrix & A, const RealTMatrix & B) { - RealTMatrix C; - C.clear(); - for (UInt i = 0; i < m; ++i) { - UInt A_i = i * k; - UInt C_i = i * n; - for (UInt j = 0; j < n; ++j) { - for (UInt l = 0; l < k; ++l) { - C[C_i + j] += A[A_i + l] * B[l * n + j]; - } - } - } - return C; -} + virtual ~Vector() { if(!wrapped) delete [] values; }; + UInt size() const { return n; }; -/* -------------------------------------------------------------------------- */ -/* vectors */ -/* -------------------------------------------------------------------------- */ + Real * storage() const { return values; }; -template -class TVector : public TMatrix { -public: - TVector(Real * data) : TMatrix(data) {}; -}; + inline Real& operator()(UInt i) { return *(values + i); }; + inline const Real& operator()(UInt i) const { return *(values + i); }; + inline Real& operator[](UInt idx) { return *(values + idx); }; + inline const Real& operator[](UInt idx) const { return *(values + idx); }; -template -class RealTVector : public TVector { -public: - RealTVector(Real * data) : TVector(data) {}; -}; + inline Vector & operator=(const Vector & vect) { + if(this != &vect) { + if(values != NULL) { + memcpy(this->values, vect.values, n * sizeof(Real)); + } else { + this->n = vect.n; + this->values = vect.values; + const_cast(vect).wrapped = true; + this->wrapped = false; + } + } + return *this; + }; -/* -------------------------------------------------------------------------- */ -/* common */ -/* -------------------------------------------------------------------------- */ + inline void clear() { memset(values, 0, n * sizeof(Real)); }; -template -inline std::ostream & operator<<(std::ostream & stream, const TMatrix & _this) -{ - _this.printself(stream); - return stream; -} + /// function to print the containt of the class + virtual void printself(std::ostream & stream, int indent = 0) const { + std::string space; + for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + stream << space << typeid(*this).name() << "<" << n <<"> :" << std::endl; + stream << space << indent << "| "; + for (UInt i = 0; i < n; ++i) { + stream << values[i] << " "; + } + stream << "|" << std::endl; + }; -/* -------------------------------------------------------------------------- */ -class Matrix { -public: - Matrix() : m(0), n(0), values(NULL), wrapped(true) {}; + // Real* &data() { return values; }; + // bool &is_wrapped() { return wrapped; }; + friend class ::akantu::Vector; - Matrix(UInt m, UInt n, Real def = 0) : m(m), n(n), values(new Real[n*m]), wrapped(false) { - std::fill_n(values, n*m, def); + protected: + UInt n; + Real* values; + bool wrapped; }; - Matrix(Real* data, UInt m, UInt n) : m(m), n(n), values(data), wrapped(true) {}; + /* -------------------------------------------------------------------------- */ + inline std::ostream & operator<<(std::ostream & stream, const Matrix & _this) + { + _this.printself(stream); + return stream; + } - Matrix(Matrix & src) { - values = src.values; - wrapped = src.wrapped; - src.wrapped = true; - }; - virtual ~Matrix() { if(!wrapped) delete [] values; }; + /* -------------------------------------------------------------------------- */ + /* templated matrix */ + /* -------------------------------------------------------------------------- */ - UInt size() { return n*m; }; + template + class TMatrix { + public: + TMatrix() : values(new T[n*m]), wrapped(false) {}; + TMatrix(T def) : values(new T[n*m]), wrapped(false) { + std::fill_n(values, n*m, def); + }; - inline Real& operator()(UInt i, UInt j) { return *(values + i*m + j); }; - inline const Real& operator()(UInt i, UInt j) const { return *(values + i*m + j); }; - inline Real& operator[](UInt idx) { return *(values + idx); }; - inline const Real& operator[](UInt idx) const { return *(values + idx); }; + TMatrix(T *data) : values(data), wrapped(true) {}; - inline Matrix & operator=(const Matrix & mat) { - if(this != &mat) { - if(values != NULL) { - memcpy(this->values, mat.values, m*n*sizeof(Real)); - } else { - this->m = mat.m; - this->n = mat.n; + TMatrix(TMatrix & src) { + values = src.values; + wrapped = src.wrapped; + src.wrapped = true; + }; - this->values = mat.values; - const_cast(mat).wrapped = true; - this->wrapped = false; - } - } - return *this; - }; + virtual ~TMatrix() { if(!wrapped) delete [] values; }; - inline Matrix operator* (const Matrix & B) { - // UInt m = this->m, n = B.n, k = this->n; - // Matrix C(m,n); - // C.clear(); + static UInt size() { return n*m; }; - // UInt A_i = 0, C_i = 0; - // for (UInt i = 0; i < m; ++i, A_i += k, C_i += n) - // for (UInt j = 0; j < n; ++j) - // for (UInt l = 0; l < k; ++l) - // C[C_i + j] += (*this)[A_i + l] * B[l * n + j]; - Matrix C(m,n); - C.mul(*this, B); - return C; - }; + inline T & operator()(UInt i, UInt j) { return *(values + i*m + j); }; + inline const T & operator()(UInt i, UInt j) const { return *(values + i*m + j); }; + inline T & operator[](UInt idx) { return *(values + idx); }; + inline const T & operator[](UInt idx) const { return *(values + idx); }; - inline void mul(const Matrix & A, const Matrix & B) { - UInt m = A.m, n = B.n, k = A.n; -#ifdef AKANTU_USE_BLAS - /// C := alpha*op(A)*op(B) + beta*C - cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, - m, n, k, - 1, - A.values, k, - B.values, n, - 0, - this->values, n); -#else - this->clear(); + inline TMatrix & operator=(const TMatrix & mat) { + if(this != &mat) { + if(values == NULL) { this->values = new T[n * m]; this->wrapped = false; } + memcpy(this->values, &mat[0], size()*sizeof(T)); + } + return *this; + }; + + inline void clear() { memset(values, 0, m * n * sizeof(T)); }; + + /// function to print the containt of the class + virtual void printself(std::ostream & stream, int indent = 0) const { + std::string space; + for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + + stream << space << typeid(*this).name() << "," << n << "," << m <<"> :" << std::endl; + for (UInt i = 0; i < m; ++i) { + stream << space << space << "| "; + for (UInt j = 0; j < n; ++j) { + stream << values[i*n +j] << " "; + } + stream << "|" << std::endl; + } + }; - UInt A_i = 0, C_i = 0; - for (UInt i = 0; i < m; ++i, A_i += k, C_i += n) - for (UInt j = 0; j < n; ++j) - for (UInt l = 0; l < k; ++l) - (*this)[C_i + j] += A[A_i + l] * B[l * n + j]; -#endif + // T* &data() { return values; }; + + friend class ::akantu::Vector; + + protected: + T * values; + bool wrapped; }; + /* -------------------------------------------------------------------------- */ - inline void clear() { memset(values, 0, m * n * sizeof(Real)); }; + template + class RealTMatrix : public TMatrix { + public: + RealTMatrix() : TMatrix() {}; + RealTMatrix(Real def) : TMatrix(def) {}; + RealTMatrix(Real * data) : TMatrix(data) {}; + RealTMatrix(RealTMatrix & mat) : TMatrix(mat) {}; + }; - /// function to print the containt of the class - virtual void printself(std::ostream & stream, int indent = 0) const { - std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - stream << space << typeid(*this).name() << "," << n << "," << m <<"> :" << std::endl; + template + inline RealTMatrix operator* (const RealTMatrix & A, const RealTMatrix & B) { + RealTMatrix C; + C.clear(); for (UInt i = 0; i < m; ++i) { - stream << space << space << "| "; + UInt A_i = i * k; + UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { - stream << values[i*n +j] << " "; + for (UInt l = 0; l < k; ++l) { + C[C_i + j] += A[A_i + l] * B[l * n + j]; + } } - stream << "|" << std::endl; } + return C; + } + + + /* -------------------------------------------------------------------------- */ + /* templated vectors */ + /* -------------------------------------------------------------------------- */ + + template + class TVector : public TMatrix { + public: + TVector(Real * data) : TMatrix(data) {}; }; - // Real* &data() { return values; }; - // bool &is_wrapped() { return wrapped; }; - friend class Vector; + template + class RealTVector : public TVector { + public: + RealTVector(Real * data) : TVector(data) {}; + }; -protected: - UInt m; - UInt n; - Real* values; - bool wrapped; -}; + /* -------------------------------------------------------------------------- */ + /* common */ + /* -------------------------------------------------------------------------- */ + template + inline std::ostream & operator<<(std::ostream & stream, const TMatrix & _this) + { + _this.printself(stream); + return stream; + } -/* -------------------------------------------------------------------------- */ -inline std::ostream & operator<<(std::ostream & stream, const Matrix & _this) -{ - _this.printself(stream); - return stream; -} +}; + __END_AKANTU__ #endif /* __AKANTU_AKA_TYPES_HH__ */ diff --git a/common/aka_types_expression.hh b/common/aka_types_expression.hh new file mode 100644 index 000000000..40b84306c --- /dev/null +++ b/common/aka_types_expression.hh @@ -0,0 +1,254 @@ +/** + * @file aka_types_expression.hh + * @author Nicolas Richart + * @date Sat Mar 26 22:20:09 2011 + * + * @brief expressions template for types operations, inspired from a work of + * Alejandro Aragón + * + * @section LICENSE + * + * Copyright (©) 2010-2011 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_types.hh" +#include "aka_math.hh" + +#ifndef __AKANTU_AKA_TYPES_EXPRESSION_HH__ +#define __AKANTU_AKA_TYPES_EXPRESSION_HH__ + +__BEGIN_AKANTU__ + +namespace types { + + /* ------------------------------------------------------------------------ */ + template struct ResultType; + + /* ------------------------------------------------------------------------ */ + template + class Expression { + public: + typedef A Type; + typedef typename A::Return Return; + typedef typename A::LeftVal LeftVal; + typedef typename A::RightVal RightVal; + + Expression(const A & a) : exp(a) { }; + + Return operator() () { return exp; }; + LeftVal left() const { return exp.left(); }; + RightVal right() const { return exp.right(); }; + + + private: + A & exp; + }; + + template + class BinaryOperation { + public: + typedef const L & LeftVal; + typedef const R & RightVal; + typedef typename ResultType > >::Type Return; + typedef Op Operator; + + BinaryOperation(LeftVal a, RightVal b) : _left(a), _right(b) {}; + BinaryOperation(L & a, RightVal b) : _left(a), _right(b) {}; + + LeftVal left() const { return _left; }; + RightVal right() const { return _right; }; + + inline Return operator() () const { + return Operator::apply(_left, _right); + } + + private: + LeftVal _left; + RightVal _right; + }; + + /* ------------------------------------------------------------------------ */ + template + class UnaryOperation { + public: + typedef const R & RightVal; + typedef const R & LeftVal; + typedef typename ResultType > >::Type Return; + typedef Op Operator; + + UnaryOperation(RightVal b) : _right(b) {}; + UnaryOperation(R & b) : _right(b) {}; + + RightVal right() const { return _right; }; + + inline Return operator() () const { + return Operator::apply(_right); + } + + private: + RightVal _right; + }; + + /* ------------------------------------------------------------------------ */ + class MultiplyOp; + class TransposeOp; + + typedef Matrix Mat; + typedef Vector Vect; + typedef Expression > MatxMat; + typedef Expression > SxMat; + typedef Expression > TrMat; + typedef Expression > TrVect; + + template <> struct ResultType > > { typedef Matrix Type; }; + + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Vector Type; }; + + template <> struct ResultType > > { typedef Matrix Type; }; + + template <> struct ResultType > > { typedef Vector Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + template <> struct ResultType > > { typedef Matrix Type; }; + + + /* ------------------------------------------------------------------------ */ + class TransposeOp { + public: + TransposeOp() {}; + }; + + /* ------------------------------------------------------------------------ */ + class MultiplyOp { + public: + MultiplyOp() {}; + + template + static inline typename ResultType > >::Type apply(const L & a, const R & b) { + return apply(a, b); + } + + + /* ---------------------------------------------------------------------- */ + /* Blas 2 */ + /* ---------------------------------------------------------------------- */ + /// @f[ y = A * x @f] + static inline Vector apply(const Mat & a, const Vect & x) { + Vector y(x.size()); + Math::matVectMul(a.rows(), a.cols(), + 1., a.storage(), x.storage(), + 0., y.storage()); + return y; + } + + /// @f[ y = A^t * x @f] + static inline Vector apply(const TrMat & at, const Vect & x) { + const Mat & a = at.right(); + Vector y(x.size()); + Math::matVectMul(a.rows(), a.cols(), + 1., a.storage(), x.storage(), + 0., y.storage()); + return y; + } + + /* ---------------------------------------------------------------------- */ + /* Blas 3 */ + /* ---------------------------------------------------------------------- */ + /// @f[ C = A * x^t @f] + static inline Mat apply(const Vect & x, const TrVect & yt) { + const Vect & y = yt.right(); + Mat c(x.size(), y.size()); + Math::matMul(x.size(), y.size(), 1, + 1., x.storage(), y.storage(), + 0., c.storage()); + return c; + } + + /// @f[ y = A^t * x^t @f] + static inline Mat apply(const TrVect & xt, const Vect & y) { + const Vect & x = xt.right(); + Mat c(y.size(), x.size()); + Math::matMul(y.size(), x.size(), 1, + 1., x.storage(), y.storage(), + 0., c.storage()); + return c; + } + + + /// @f[ C = A * B @f] + static inline Mat apply(const Mat & a, const Mat & b) { + Mat c(a.rows(), b.cols()); + Math::matMul(a.rows(), b.cols(), a.cols(), + 1., a.storage(), b.storage(), + 0., c.storage()); + return c; + } + + /// @f[ C = A^t * B @f] + static inline Mat apply(const TrMat & at, const Mat & b) { + const Mat & a = at.right(); + Mat c(a.rows(), b.cols()); + Math::matMul(a.rows(), b.cols(), a.cols(), + 1., a.storage(), b.storage(), + 0., c.storage()); + return c; + } + + /// @f[ C = A * B^t @f] + static inline Mat apply(const Mat & a, const TrMat & bt) { + const Mat & b = bt.right(); + Mat c(a.rows(), b.cols()); + Math::matMul(a.rows(), b.cols(), a.cols(), + 1., a.storage(), b.storage(), + 0., c.storage()); + return c; + } + + /// @f[ C = A^t * B^t @f] + static inline Mat apply(const TrMat & at, const TrMat & bt) { + const Mat & a = at.right(); + const Mat & b = bt.right(); + Mat c(a.rows(), b.cols()); + Math::matMul(a.rows(), b.cols(), a.cols(), + 1., a.storage(), b.storage(), + 0., c.storage()); + return c; + } + }; + + template + typename ResultType > >::Type operator* (const A & a, const A & b) { + return Expression >(BinaryOperation(a,b)); + } + +}; // namespace types + +__END_AKANTU__ + +#endif /* __AKANTU_AKA_TYPES_EXPRESSION_HH__ */ diff --git a/common/aka_vector.hh b/common/aka_vector.hh index 109d785a2..301b8e77f 100644 --- a/common/aka_vector.hh +++ b/common/aka_vector.hh @@ -1,268 +1,270 @@ /** * @file aka_vector.hh * @author Nicolas Richart * @date Thu Jun 17 10:04:55 2010 * * @brief class of vectors * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_VECTOR_HH__ #define __AKANTU_VECTOR_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Matrix; /// class that afford to store vectors in static memory class VectorBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: VectorBase(const VectorID & id = ""); virtual ~VectorBase(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the amount of space allocated in bytes inline UInt getMemorySize() const; /// set the size to zero without freeing the allocated space inline void empty(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt); AKANTU_GET_MACRO(Size, size, UInt); AKANTU_GET_MACRO(NbComponent, nb_component, UInt); AKANTU_GET_MACRO(ID, id, const VectorID &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the vector VectorID id; /// the size allocated UInt allocated_size; /// the size used UInt size; /// number of components UInt nb_component; /// size of the stored type UInt size_of_type; }; - +namespace types { + class Matrix; +}; /* -------------------------------------------------------------------------- */ template class Vector : public VectorBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef T value_type; typedef value_type & reference; typedef value_type * pointer_type; typedef const value_type & const_reference; /// Allocation of a new vector Vector(UInt size = 0, UInt nb_component = 1, const VectorID & id = ""); /// Allocation of a new vector with a default value Vector(UInt size, UInt nb_component, const value_type def_values[], const VectorID & id = ""); /// Allocation of a new vector with a default value Vector(UInt size, UInt nb_component, const_reference value, const VectorID & id = ""); /// Copy constructor (deep copy if deep=true) \todo to implement Vector(const Vector& vect, bool deep = true); virtual ~Vector(); /* ------------------------------------------------------------------------ */ /* Iterator */ /* ------------------------------------------------------------------------ */ public: template class iterator { public: typedef R returned_type; typedef returned_type & ref_returned_type; iterator(); iterator(pointer_type data, UInt offset); iterator(const iterator & it); ~iterator(); inline iterator & operator=(const iterator & it); inline ref_returned_type operator*() { return *ret; }; inline iterator & operator++(); inline iterator & operator+=(const UInt n); inline bool operator==(const iterator & other); inline bool operator!=(const iterator & other); private: UInt offset; returned_type * ret; }; template inline iterator begin(); template inline iterator end(); - inline iterator begin(UInt m, UInt n); - inline iterator end(UInt m, UInt n); + inline iterator begin(UInt m, UInt n); + inline iterator end(UInt m, UInt n); inline reference operator()(UInt i, UInt j = 0); inline const_reference operator()(UInt i, UInt j = 0) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get jth componemt of the ith tuple in read-only inline const_reference get(UInt i, UInt j = 0) const; /// get jth componemt of the ith tuple inline reference at(UInt i, UInt j = 0); /// add an element at the and of the vector with the value value for all /// component inline void push_back(const_reference value); /// add an element at the and of the vector inline void push_back(const value_type new_elem[]); /** * remove an element and move the last one in the hole * /!\ change the order in the vector */ inline void erase(UInt i); /// change the size of the vector and allocate more memory if needed void resize(UInt size); /// change the number of components by interlacing data void extendComponentsInterlaced(UInt multiplicator, UInt stride); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; // Vector& operator=(const Vector& vect); /// search elem in the vector, return the position of the first occurrence or /// -1 if not found Int find(const_reference elem) const; /// set a vvector to 0 inline void clear() { memset(values, 0, size*nb_component*sizeof(T)); }; /// copy the content of an other vector void copy(const Vector & vect); protected: /// perform the allocation for the constructors void allocate(UInt size, UInt nb_component = 1); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: UInt getSize() const{ return this->size; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: /// array of values T *values; // /!\ very dangerous }; __END_AKANTU__ #include "aka_types.hh" __BEGIN_AKANTU__ #include "aka_vector_inline_impl.cc" /* -------------------------------------------------------------------------- */ /* Inline Functions Vector */ /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Vector & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ /* Inline Functions VectorBase */ /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const VectorBase & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_VECTOR_HH__ */ diff --git a/common/aka_vector_inline_impl.cc b/common/aka_vector_inline_impl.cc index e29ba8be6..5dca50931 100644 --- a/common/aka_vector_inline_impl.cc +++ b/common/aka_vector_inline_impl.cc @@ -1,279 +1,279 @@ /** * @file aka_vector_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 15 00:09:33 2010 * * @brief Inline functions of the classes Vector and VectorBase * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Vector */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline T & Vector::operator()(UInt i, UInt j) { AKANTU_DEBUG_ASSERT(size > 0, "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range"); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline const T & Vector::operator()(UInt i, UInt j) const { AKANTU_DEBUG_ASSERT(size > 0, "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range"); return values[i*nb_component + j]; } template inline T & Vector::at(UInt i, UInt j) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(size > 0, "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range"); AKANTU_DEBUG_OUT(); return values[i*nb_component + j]; } template inline const T & Vector::get(UInt i, UInt j) const{ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(size > 0, "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range"); AKANTU_DEBUG_OUT(); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline void Vector::push_back(const T & value) { // AKANTU_DEBUG_IN(); UInt pos = size; resize(size+1); /// @todo see if with std::uninitialized_fill it allow to build vector of objects for (UInt i = 0; i < nb_component; ++i) { values[pos*nb_component + i] = value; } // AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Vector::push_back(const T new_elem[]) { // AKANTU_DEBUG_IN(); UInt pos = size; resize(size+1); for (UInt i = 0; i < nb_component; ++i) { values[pos*nb_component + i] = new_elem[i]; } // AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Vector::erase(UInt i){ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((size > 0), "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size), "The element at position [" << i << "] is out of range"); if(i != (size - 1)) { for (UInt j = 0; j < nb_component; ++j) { values[i*nb_component + j] = values[(size-1)*nb_component + j]; } } resize(size - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline Vector::iterator Vector::begin() { return iterator(values, nb_component); } /* -------------------------------------------------------------------------- */ template template inline Vector::iterator Vector::end() { return iterator(values + nb_component * size, nb_component); } /* -------------------------------------------------------------------------- */ /* Inline Functions VectorBase */ /* -------------------------------------------------------------------------- */ inline UInt VectorBase::getMemorySize() const { return allocated_size * nb_component * size_of_type; } inline void VectorBase::empty() { size = 0; } /* -------------------------------------------------------------------------- */ /* Iterators */ /* -------------------------------------------------------------------------- */ template template Vector::iterator::iterator() : offset(0), ret(NULL) { } /* -------------------------------------------------------------------------- */ template template Vector::iterator::iterator(pointer_type data, UInt offset) : offset(offset), ret(new returned_type(data)) { AKANTU_DEBUG_ASSERT(offset == ret->size(), "The iterator is not compatible with the type " << typeid(returned_type).name()); } /* -------------------------------------------------------------------------- */ template template Vector::iterator::iterator(const Vector::iterator & it) { if(this != &it) { this->offset = it.offset; this->ret = new returned_type(it.ret->values); } } /* -------------------------------------------------------------------------- */ template template Vector::iterator::~iterator() { delete ret; } /* -------------------------------------------------------------------------- */ template template inline Vector::iterator & Vector::iterator::operator++() { ret->values += offset; return *this; } /* -------------------------------------------------------------------------- */ template template inline Vector::iterator & Vector::iterator::operator=(const Vector::iterator & it) { if(this != &it) { this->offset = it.offset; this->ret = new returned_type(it.ret->values); } return *this; } /* -------------------------------------------------------------------------- */ template template inline Vector::iterator & Vector::iterator::operator+=(const UInt n) { ret->values += offset * n; return *this; } /* -------------------------------------------------------------------------- */ template template inline bool Vector::iterator::operator==(const iterator & other) { return (*this).ret->values == other.ret->values; } /* -------------------------------------------------------------------------- */ template template inline bool Vector::iterator::operator!=(const iterator & other) { return (*this).ret->values != other.ret->values; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template<> template<> -class Vector::iterator { +class Vector::iterator { public: iterator() : ret(NULL), offset(0) { }; - iterator(Real * data, UInt offset, UInt m, UInt n) : ret(new Matrix(data, m, n)), offset(offset) { + iterator(Real * data, UInt offset, UInt m, UInt n) : ret(new types::Matrix(data, m, n)), offset(offset) { AKANTU_DEBUG_ASSERT(offset == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); }; iterator(const iterator & it) { if(this != &it) { offset = it.offset; - ret = new Matrix(it.ret->values, it.ret->m, it.ret->n); + ret = new types::Matrix(it.ret->values, it.ret->m, it.ret->n); } }; ~iterator() { delete ret; }; inline iterator & operator=(const iterator & it) { if(this != &it) { offset = it.offset; - ret = new Matrix(it.ret->values, it.ret->m, it.ret->n); + ret = new types::Matrix(it.ret->values, it.ret->m, it.ret->n); } return *this; }; - inline Matrix & operator*() { return *ret; }; + inline types::Matrix & operator*() { return *ret; }; inline iterator & operator++() { ret->values += offset; return *this; }; inline iterator & operator+=(const UInt n) { ret->values += n*offset; return *this; }; inline bool operator==(const iterator & other) { return ret->values == other.ret->values; }; inline bool operator!=(const iterator & other) { return ret->values != other.ret->values; }; private: - Matrix * ret; - UInt offset; + types::Matrix * ret; + UInt offset; }; /* -------------------------------------------------------------------------- */ template<> -inline Vector::iterator Vector::begin(UInt m, UInt n) { - return iterator(values, nb_component, m, n); +inline Vector::iterator Vector::begin(UInt m, UInt n) { + return iterator(values, nb_component, m, n); } /* -------------------------------------------------------------------------- */ template<> -inline Vector::iterator Vector::end(UInt m, UInt n) { - return iterator(values + nb_component * size, nb_component, m, n); +inline Vector::iterator Vector::end(UInt m, UInt n) { + return iterator(values + nb_component * size, nb_component, m, n); } diff --git a/fem/fem.hh b/fem/fem.hh index 2efcdf9ae..cd65d9ff8 100644 --- a/fem/fem.hh +++ b/fem/fem.hh @@ -1,221 +1,223 @@ /** * @file fem.hh * @author Nicolas Richart * @date Fri Jul 16 10:24:24 2010 * * @brief FEM class * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_FEM_HH__ #define __AKANTU_FEM_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "mesh.hh" #include "element_class.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class FEM : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FEM(Mesh & mesh, UInt spatial_dimension = 0, FEMID id = "fem", MemoryID memory_id = 0); virtual ~FEM(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile of the sparse matrix corresponding to the mesh void initSparseMatrixProfile(SparseMatrixType sparse_matrix_type = _unsymmetric); /// pre-compute all the shape functions, their derivatives and the jacobians virtual void initShapeFunctions(GhostType ghost_type = _not_ghost)=0; /* ------------------------------------------------------------------------ */ /* Integration method bridges */ /* ------------------------------------------------------------------------ */ /// integrate f for all elements of type "type" virtual void integrate(const Vector & f, Vector &intf, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const = 0; /// integrate a scalar value on all elements of type "type" virtual Real integrate(const Vector & f, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const = 0; /* ------------------------------------------------------------------------ */ /* compatibility with old FEM fashion */ /* ------------------------------------------------------------------------ */ /// get the number of quadrature points virtual UInt getNbQuadraturePoints(const ElementType & type)=0; /// get the precomputed shapes const virtual Vector & getShapes(const ElementType & type)=0; /// get the precomputed ghost shapes const virtual Vector & getGhostShapes(const ElementType & type)=0; /// get the derivatives of shapes const virtual Vector & getShapesDerivatives(const ElementType & type)=0; /// get the ghost derivatives of shapes const virtual Vector & getGhostShapesDerivatives(const ElementType & type)=0; /// get quadrature points const virtual Vector & getQuadraturePoints(const ElementType & type)=0; /* ------------------------------------------------------------------------ */ /* Shape method bridges */ /* ------------------------------------------------------------------------ */ virtual void gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, const UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL)=0; virtual void interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const =0; /* ------------------------------------------------------------------------ */ /* Other methods */ /* ------------------------------------------------------------------------ */ /// pre-compute normals on control points virtual void computeNormalsOnControlPoints(GhostType ghost_type = _not_ghost)=0; /// assemble vectors void assembleVector(const Vector & elementary_vect, Vector & nodal_values, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL, Real scale_factor = 1) const; /// assemble matrix in the complete sparse matrix void assembleMatrix(const Vector & elementary_mat, SparseMatrix & matrix, const Vector & equation_number, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; - virtual void assembleFieldLumped(const ByElementTypeReal & field_1, + /// assemble a field as a lumped matrix (ex. rho in lumped mass) + virtual void assembleFieldLumped(const Vector & field_1, Vector & lumped, + ElementType type, GhostType ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; private: /// initialise the class void init(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt); /// get the mesh contained in the fem object inline Mesh & getMesh() const; /// get the in-radius of an element static inline Real getElementInradius(Real * coord, const ElementType & type); /// get the normals on quadrature points AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NormalsOnQuadPoints, normals_on_quad_points, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the fem object FEMID id; /// spatial dimension of the problem UInt element_dimension; /// the mesh on which all computation are made Mesh * mesh; /// normals at quadrature points ByElementTypeReal normals_on_quad_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "fem_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const FEM & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "fem_template.hh" #endif /* __AKANTU_FEM_HH__ */ diff --git a/fem/fem_template.cc b/fem/fem_template.cc index b48e2a0d2..103cc7ac2 100644 --- a/fem/fem_template.cc +++ b/fem/fem_template.cc @@ -1,518 +1,513 @@ /** * @file fem_template.cc * @author Guillaume Anciaux * @author Nicolas Richart * @date Fri Feb 11 11:37:47 2011 * * @brief implementation of the generic FEMTemplate class * * @section LICENSE * * Copyright (©) 2010-2011 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 "integrator_gauss.hh" #include "shape_lagrange.hh" #include "fem.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template FEMTemplate::FEMTemplate(Mesh & mesh, UInt spatial_dimension, FEMID id,MemoryID memory_id) :FEM(mesh,spatial_dimension,id,memory_id), integrator(mesh), shape_functions(mesh) { } /* -------------------------------------------------------------------------- */ template FEMTemplate::~FEMTemplate() { } /* -------------------------------------------------------------------------- */ template void FEMTemplate::gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, const UInt nb_degree_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements){ AKANTU_DEBUG_IN(); #define COMPUTE_GRADIENT(type) \ if (element_dimension == ElementClass::getSpatialDimension()) \ shape_functions.template gradientOnControlPoints(u, \ nablauq, \ nb_degree_of_freedom, \ ghost_type, \ filter_elements); AKANTU_BOOST_ELEMENT_SWITCH(COMPUTE_GRADIENT) #undef COMPUTE_GRADIENT AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FEMTemplate::initShapeFunctions(GhostType ghost_type){ AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh->getSpatialDimension(); const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; #define INIT_SHAPE_FUNCTIONS(type) \ if (element_dimension != ElementClass::getSpatialDimension()) continue; \ integrator.template computeQuadraturePoints(); \ integrator.template precomputeJacobiansOnQuadraturePoints(spatial_dimension,ghost_type); \ Vector & control_points = integrator.template getQuadraturePoints(); \ shape_functions.template setControlPointsByType(control_points); \ shape_functions.template precomputeShapesOnControlPoints(ghost_type); \ if (element_dimension == spatial_dimension) \ shape_functions.template precomputeShapeDerivativesOnControlPoints(ghost_type); AKANTU_BOOST_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS) #undef INIT_SHAPE_FUNCTIONS } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FEMTemplate::integrate(const Vector & f, Vector &intf, UInt nb_degree_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const{ #define INTEGRATE(type) \ integrator.template integrate(f, \ intf, \ nb_degree_of_freedom, \ ghost_type, \ filter_elements); AKANTU_BOOST_ELEMENT_SWITCH(INTEGRATE) #undef INTEGRATE } /* -------------------------------------------------------------------------- */ template Real FEMTemplate::integrate(const Vector & f, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const{ AKANTU_DEBUG_IN(); Real integral = 0.; #define INTEGRATE(type) \ integral = integrator.template integrate(f, \ ghost_type, \ filter_elements); AKANTU_BOOST_ELEMENT_SWITCH(INTEGRATE) #undef INTEGRATE AKANTU_DEBUG_OUT(); return integral; } /* -------------------------------------------------------------------------- */ template void FEMTemplate::interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degree_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const{ AKANTU_DEBUG_IN(); #define INTERPOLATE(type) \ shape_functions.template interpolateOnControlPoints(u, \ uq, \ nb_degree_of_freedom, \ ghost_type, \ filter_elements); AKANTU_BOOST_ELEMENT_SWITCH(INTERPOLATE) #undef INTERPOLATE AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FEMTemplate::computeNormalsOnControlPoints(GhostType ghost_type) { AKANTU_DEBUG_IN(); // Real * coord = mesh->getNodes().values; UInt spatial_dimension = mesh->getSpatialDimension(); //allocate the normal arrays if (ghost_type == _not_ghost) mesh->initByElementTypeRealVector(normals_on_quad_points,spatial_dimension,element_dimension, id,"normals_onquad",ghost_type); else{ AKANTU_DEBUG_ERROR("to be implemented"); } //loop over the type to build the normals const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; UInt element_type_spatial_dimension = Mesh::getSpatialDimension(type); if(element_type_spatial_dimension != element_dimension) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quad_points = getQuadraturePoints(type).getSize(); UInt * elem_val; UInt nb_element; Real * normals_on_quad_val = NULL; if(ghost_type == _not_ghost) { elem_val = mesh->getConnectivity(type).values; nb_element = mesh->getConnectivity(type).getSize(); normals_on_quad_points[type]->resize(nb_element * nb_quad_points); normals_on_quad_val = normals_on_quad_points[type]->values; } else { elem_val = mesh->getGhostConnectivity(type).values; nb_element = mesh->getGhostConnectivity(type).getSize(); } /* ---------------------------------------------------------------------- */ #define COMPUTE_NORMALS_ON_QUAD(type) \ do { \ Vector & quads = integrator. template getQuadraturePoints(); \ UInt nb_points = quads.getSize(); \ Real local_coord[spatial_dimension * nb_nodes_per_element]; \ for (UInt elem = 0; elem < nb_element; ++elem) { \ mesh->extractNodalCoordinatesFromElement(local_coord, \ elem_val+elem*nb_nodes_per_element, \ nb_nodes_per_element); \ ElementClass::computeNormalsOnQuadPoint(local_coord, \ spatial_dimension, \ normals_on_quad_val); \ normals_on_quad_val += spatial_dimension*nb_points; \ } \ } while(0) /* ---------------------------------------------------------------------- */ AKANTU_BOOST_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_QUAD) #undef COMPUTE_NORMALS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* compatibility functions */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline UInt FEMTemplate::getNbQuadraturePoints(const ElementType & type) { AKANTU_DEBUG_IN(); UInt nb_quad_points = 0; #define GET_NB_QUAD(type) \ nb_quad_points = integrator. template getQuadraturePoints().getSize(); AKANTU_BOOST_ELEMENT_SWITCH(GET_NB_QUAD) #undef GET_NB_QUAD AKANTU_DEBUG_OUT(); return nb_quad_points; } /* -------------------------------------------------------------------------- */ template inline const Vector & FEMTemplate::getShapes(const ElementType & type) { AKANTU_DEBUG_IN(); const Vector * ret = NULL; #define GET_SHAPES(type) \ ret = &(shape_functions.getShapes(type)); AKANTU_BOOST_ELEMENT_SWITCH(GET_SHAPES) #undef GET_SHAPES AKANTU_DEBUG_OUT(); return *ret; } /* -------------------------------------------------------------------------- */ template inline const Vector & FEMTemplate::getGhostShapes(const ElementType & type) { AKANTU_DEBUG_IN(); const Vector * ret = NULL; #define GET_SHAPES(type) \ ret = &(shape_functions.getGhostShapes(type)); AKANTU_BOOST_ELEMENT_SWITCH(GET_SHAPES) #undef GET_SHAPES AKANTU_DEBUG_OUT(); return *ret; } /* -------------------------------------------------------------------------- */ template inline const Vector & FEMTemplate::getShapesDerivatives(const ElementType & type) { AKANTU_DEBUG_IN(); const Vector * ret = NULL; #define GET_SHAPES(type) \ ret = &(shape_functions.getShapesDerivatives(type)); AKANTU_BOOST_ELEMENT_SWITCH(GET_SHAPES) #undef GET_SHAPES AKANTU_DEBUG_OUT(); return *ret; } /* -------------------------------------------------------------------------- */ template inline const Vector & FEMTemplate::getGhostShapesDerivatives(const ElementType & type) { AKANTU_DEBUG_IN(); const Vector * ret = NULL; #define GET_SHAPES(type) \ ret = &(shape_functions.getGhostShapesDerivatives(type)); AKANTU_BOOST_ELEMENT_SWITCH(GET_SHAPES) #undef GET_SHAPES AKANTU_DEBUG_OUT(); return *ret; } /* -------------------------------------------------------------------------- */ template inline const Vector & FEMTemplate::getQuadraturePoints(const ElementType & type) { AKANTU_DEBUG_IN(); const Vector * ret = NULL; #define GET_QUADS(type) \ ret = &(integrator. template getQuadraturePoints()); AKANTU_BOOST_ELEMENT_SWITCH(GET_QUADS) #undef GET_QUADS AKANTU_DEBUG_OUT(); return *ret; } /* -------------------------------------------------------------------------- */ /* Matrix lumping functions */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template -void FEMTemplate::assembleFieldLumped(const ByElementTypeReal & field_1, +void FEMTemplate::assembleFieldLumped(const Vector & field_1, Vector & lumped, + ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(ghost_type); - Mesh::ConnectivityTypeList::const_iterator it; #define ASSEMBLE_LUMPED(type) \ - assembleLumpedTemplate(*field_1[type], lumped, ghost_type) + assembleLumpedTemplate(field_1, lumped, ghost_type) - for(it = type_list.begin(); it != type_list.end(); ++it) { - if(Mesh::getSpatialDimension(*it) != element_dimension) continue; - ElementType type = *it; - AKANTU_BOOST_ELEMENT_SWITCH(ASSEMBLE_LUMPED); - } + AKANTU_BOOST_ELEMENT_SWITCH(ASSEMBLE_LUMPED); #undef ASSEMBLE_LUMPED AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void FEMTemplate::assembleLumpedTemplate(const Vector & field_1, Vector & lumped, GhostType ghost_type) { this->template assembleLumpedRowSum(field_1, lumped, ghost_type); } /* -------------------------------------------------------------------------- */ //template template <> template <> void FEMTemplate::assembleLumpedTemplate<_triangle_6>(const Vector & field_1, Vector & lumped, GhostType ghost_type) { assembleLumpedDiagonalScaling<_triangle_6>(field_1, lumped, ghost_type); } // /* -------------------------------------------------------------------------- */ template <> template <> void FEMTemplate::assembleLumpedTemplate<_tetrahedron_10>(const Vector & field_1, Vector & lumped, GhostType ghost_type) { assembleLumpedDiagonalScaling<_tetrahedron_10>(field_1, lumped, ghost_type); } /* -------------------------------------------------------------------------- */ /** * @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$ */ template template void FEMTemplate::assembleLumpedRowSum(const Vector & field_1, Vector & lumped, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt shapes_size = ElementClass::getShapeSize(); Vector * field_times_shapes = new Vector(0, 1);//shapes_size); shape_functions.template fieldTimesShapes(field_1, *field_times_shapes, ghost_type); Vector * int_field_times_shapes = new Vector(0, shapes_size, "inte_rho_x_shapes"); integrator.template integrate(*field_times_shapes, *int_field_times_shapes, shapes_size, ghost_type, NULL); delete field_times_shapes; assembleVector(*int_field_times_shapes, lumped, 1, type, ghost_type); delete int_field_times_shapes; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$ */ template template void FEMTemplate::assembleLumpedDiagonalScaling(const Vector & field_1, Vector & lumped, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element_p1 = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type)); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = integrator.template getQuadraturePoints().getSize(); UInt nb_element = field_1.getSize() / nb_quadrature_points; Real corner_factor = 0; Real mid_factor = 0; switch(type){ case _triangle_6 : corner_factor = 1./12.; mid_factor = 1./4.; break; case _tetrahedron_10: corner_factor = 1./32.; mid_factor = 7./48.; break; default: corner_factor = 0; mid_factor = 0; } if (nb_element == 0) { AKANTU_DEBUG_OUT(); return; } /// compute @f$ \int \rho dV = \rho V @f$ for each element Vector * int_field_1 = new Vector(field_1.getSize(), 1, "inte_rho_x_1"); integrator.template integrate(field_1, *int_field_1, 1, ghost_type, NULL); /// distribute the mass of the element to the nodes Vector * lumped_per_node = new Vector(nb_element, nb_nodes_per_element, "mass_per_node"); Real * int_field_1_val = int_field_1->values; Real * lumped_per_node_val = lumped_per_node->values; for (UInt e = 0; e < nb_element; ++e) { Real lmass = *int_field_1_val * corner_factor; for (UInt n = 0; n < nb_nodes_per_element_p1; ++n) *lumped_per_node_val++ = lmass; /// corner points lmass = *int_field_1_val * mid_factor; for (UInt n = nb_nodes_per_element_p1; n < nb_nodes_per_element; ++n) *lumped_per_node_val++ = lmass; /// mid points int_field_1_val++; } delete int_field_1; assembleVector(*lumped_per_node, lumped, 1, type, ghost_type); delete lumped_per_node; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* template instanciation */ /* -------------------------------------------------------------------------- */ template class FEMTemplate; __END_AKANTU__ diff --git a/fem/fem_template.hh b/fem/fem_template.hh index 2f3353c39..f3118b9b0 100644 --- a/fem/fem_template.hh +++ b/fem/fem_template.hh @@ -1,176 +1,177 @@ /** * @file fem_template.hh * @author Guillaume Anciaux * @date Thu Feb 10 10:55:21 2011 * * @brief templated class that calls integration and shape objects * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ #ifndef __AKANTU_FEM_TEMPLATE_HH__ #define __AKANTU_FEM_TEMPLATE_HH__ /* -------------------------------------------------------------------------- */ #include "fem.hh" #include "integrator.hh" #include "shape_functions.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class FEMTemplate : public FEM{ /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FEMTemplate(Mesh & mesh, UInt spatial_dimension = 0, FEMID id = "fem", MemoryID memory_id = 0); virtual ~FEMTemplate(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// pre-compute all the shape functions, their derivatives and the jacobians void initShapeFunctions(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Integration method bridges */ /* ------------------------------------------------------------------------ */ /// integrate f for all elements of type "type" void integrate(const Vector & f, Vector &intf, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// integrate a scalar value on all elements of type "type" Real integrate(const Vector & f, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// get the number of quadrature points UInt getNbQuadraturePoints(const ElementType & type); /// get shapes precomputed const Vector & getShapes(const ElementType & type); /// get the ghost shapes precomputed const Vector & getGhostShapes(const ElementType & type); /// get the derivatives of shapes const Vector & getShapesDerivatives(const ElementType & type); /// get the ghost derivatives of shapes const Vector & getGhostShapesDerivatives(const ElementType & type); /// get quadrature points const Vector & getQuadraturePoints(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Shape method bridges */ /* ------------------------------------------------------------------------ */ void gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, const UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL); void interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /* ------------------------------------------------------------------------ */ /* Other methods */ /* ------------------------------------------------------------------------ */ /// pre-compute normals on control points void computeNormalsOnControlPoints(GhostType ghost_type = _not_ghost); /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const{}; - void assembleFieldLumped(const ByElementTypeReal & field_1, + void assembleFieldLumped(const Vector & field_1, Vector & lumped, + ElementType type, GhostType ghost_type); private: template void assembleLumpedTemplate(const Vector & field_1, Vector & lumped, GhostType ghost_type); template void assembleLumpedRowSum(const Vector & field_1, Vector & lumped, GhostType ghost_type); template void assembleLumpedDiagonalScaling(const Vector & field_1, Vector & lumped, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: Integ integrator; Shape shape_functions; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "fem_template_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const FEMTemplate & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_FEM_TEMPLATE_HH__ */ diff --git a/model/solid_mechanics/solid_mechanics_model.cc b/model/solid_mechanics/solid_mechanics_model.cc index 8288a0766..1bb28e52e 100644 --- a/model/solid_mechanics/solid_mechanics_model.cc +++ b/model/solid_mechanics/solid_mechanics_model.cc @@ -1,727 +1,731 @@ /** * @file solid_mechanics_model.cc * @author Nicolas Richart * @date Thu Jul 22 14:35:38 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 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 "solid_mechanics_model.hh" #include "aka_math.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "sparse_matrix.hh" #include "solver.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ModelID & id, const MemoryID & memory_id) : Model(id, memory_id), time_step(NAN), f_m2a(1.0), + mass_matrix(NULL), + velocity_damping_matrix(NULL), stiffness_matrix(NULL), integrator(new CentralDifference()), increment_flag(false), solver(NULL), spatial_dimension(dim) { AKANTU_DEBUG_IN(); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; for(UInt t = _not_defined; t < _max_element_type; ++t) { this->element_material[t] = NULL; this->ghost_element_material[t] = NULL; } registerTag(_gst_smm_mass, "Mass"); registerTag(_gst_smm_for_strain, "Explicit Residual"); registerTag(_gst_smm_boundary, "Boundary conditions"); materials.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; if(solver) delete solver; + if(mass_matrix) delete mass_matrix; + if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix) delete stiffness_matrix; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); mass = &(alloc(sstr_mass.str(), nb_nodes, 1)); // \todo see if it must not be spatial_dimension velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = getFEM().getMesh().getNbElement(*it); if(!element_material[*it]) { std::stringstream sstr_elma; sstr_elma << id << ":element_material:" << *it; element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } } const Mesh::ConnectivityTypeList & ghost_type_list = getFEM().getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { UInt nb_element = getFEM().getMesh().getNbGhostElement(*it); std::stringstream sstr_elma; sstr_elma << id << ":ghost_element_material:" << *it; ghost_element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); current_position->resize(nb_nodes); //Vector * current_position = new Vector(nb_nodes, spatial_dimension, NAN, "position"); Real * current_position_val = current_position->values; Real * position_val = getFEM().getMesh().getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); if (need_initialize) initializeUpdateResidualData(); /// start synchronization asynchronousSynchronize(_gst_smm_for_strain); /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { //(*mat_it)->updateResidual(*current_position, _not_ghost); (*mat_it)->updateResidual(*displacement, _not_ghost); } /// finalize communications waitEndSynchronize(_gst_smm_for_strain); /// call update residual on each ghost elements for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { //(*mat_it)->updateResidual(*current_position, _ghost); (*mat_it)->updateResidual(*displacement, _ghost); } // current_position; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); UInt nb_nodes = acceleration->getSize(); UInt nb_degre_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * residual_val = residual->values; Real * accel_val = acceleration->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { for (UInt d = 0; d < nb_degre_of_freedom; d++) { if(!(*boundary_val)) { *accel_val = f_m2a * *residual_val / *mass_val; } residual_val++; accel_val++; boundary_val++; } mass_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real)); } integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_nodes = displacement->getSize(); for (UInt n = 0; n < nb_nodes; ++n) { *inc_val = *dis_val - *inc_val; *inc_val++; *dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorr(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initImplicitSolver() { AKANTU_DEBUG_IN(); const Mesh & mesh = getFEM().getMesh(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); std::stringstream sstr_eq; sstr_eq << id << ":equation_number"; equation_number = &(alloc(sstr_eq.str(), nb_nodes, spatial_dimension, 0)); Int * equation_number_val = equation_number->values; for (UInt n = 0; n < nb_nodes; ++n) { UInt real_n = mesh.getNodeGlobalId(n); UInt is_local_node = mesh.isLocalNode(n); for (UInt d = 0; d < spatial_dimension; ++d) { UInt global_eq_num = (is_local_node ? real_n : nb_global_node + real_n) * spatial_dimension + d; *(equation_number_val++) = global_eq_num; local_eq_num_to_global[global_eq_num] = n * spatial_dimension + d; } } stiffness_matrix->buildProfile(mesh, *equation_number); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver_stiffness_matrix"; solver = new SolverMumps(*stiffness_matrix, sstr_solv.str()); dynamic_cast(solver)->initNodesLocation(getFEM().getMesh(), spatial_dimension); solver->initialize(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); updateCurrentPosition(); /// start synchronization // asynchronousSynchronize(_gst_smm_for_strain); stiffness_matrix->clear(); /// call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(*current_position, _not_ghost); } /// finalize communications // waitEndSynchronize(_gst_smm_for_strain); // /// call compute stiffness matrix on each ghost elements // for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { // (*mat_it)->computeStiffnessMatrix(*current_position, _ghost); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solve() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving an implicit step."); UInt nb_nodes = displacement->getSize(); UInt nb_degre_of_freedom = displacement->getNbComponent(); stiffness_matrix->applyBoundary(*boundary, local_eq_num_to_global); // stiffness_matrix->removeBoundary(*boundary); // #ifdef AKANTU_USE_MUMPS // std::stringstream sstr_solv; sstr_solv << id << ":solver_stiffness_matrix"; // solver = new SolverMumps(*stiffness_matrix, sstr_solv.str()); // dynamic_cast(solver)->initNodesLocation(getFEM().getMesh(), nb_degre_of_freedom); // #else // AKANTU_DEBUG_ERROR("You should at least activate one solver."); // #endif //AKANTU_USE_MUMPS // solver->initialize(); // Vector * tmp = new Vector(0, 1); // for (UInt i = 0; i < nb_degre_of_freedom * nb_nodes; ++i) { // if(! boundary->values[i]) { // tmp->push_back(residual->values[i]); // } // } // solver->setRHS(*tmp); // tmp->clear(); // not needed in fact // solver->solve(*tmp); solver->setRHS(*residual); if(!increment) setIncrementFlagOn(); solver->solve(*increment); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; // Real * tmp_val = tmp->values; for (UInt n = 0; n < nb_nodes * nb_degre_of_freedom; ++n) { if(!(*boundary_val)) { // *increment_val = *(tmp_val++); *displacement_val += *increment_val; } displacement_val++; boundary_val++; increment_val++; } //stiffness_matrix->restoreProfile(); // delete tmp; // delete solver; // solver = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) { AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degre_of_freedom = displacement->getNbComponent(); Real norm = 0; Real * increment_val = increment->values; bool * boundary_val = boundary->values; Mesh & mesh = getFEM().getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degre_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm += *increment_val * *increment_val; } boundary_val++; increment_val++; } } allReduce(&norm, _so_sum); AKANTU_DEBUG_INFO("Norm of increment : " << sqrt(norm)); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (sqrt(norm) < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); Real norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; Mesh & mesh = getFEM().getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } allReduce(&norm, _so_sum); AKANTU_DEBUG_INFO("Norm of residual : " << sqrt(norm)); AKANTU_DEBUG_ASSERT(!isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (sqrt(norm) < tolerance); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); Real * coord = getFEM().getMesh().getNodes().values; Real * disp_val = displacement->values; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = getFEM().getMesh().getNbNodesPerElement(*it); UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt * conn = getFEM().getMesh().getConnectivity(*it).values; UInt * elem_mat_val = element_material[*it]->values; Real * u = new Real[nb_nodes_per_element*spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(u + n * spatial_dimension, coord + offset_conn, spatial_dimension * sizeof(Real)); for (UInt i = 0; i < spatial_dimension; ++i) { u[n * spatial_dimension + i] += disp_val[offset_conn + i]; } } Real el_size = getFEM().getElementInradius(u, *it); Real el_dt = mat_val[elem_mat_val[el]]->getStableTimeStep(el_size); min_dt = min_dt > el_dt ? el_dt : min_dt; } delete [] u; } /// reduction min over all processors allReduce(&min_dt, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { epot += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors allReduce(&epot, _so_sum); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = getFEM().getMesh().getNbNodes(); // Vector * v_square = new Vector(nb_nodes, 1, "v_square"); Real * vel_val = velocity->values; Real * mass_val = mass->values; Mesh & mesh = getFEM().getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { Real v2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt i = 0; i < spatial_dimension; ++i) { if(is_local_node) { v2 += *vel_val * *vel_val; } vel_val++; } ekin += *mass_val * v2; mass_val++; } // Real * v_s_val = v_square->values; // for (UInt n = 0; n < nb_nodes; ++n) { // *v_s_val = 0; // for (UInt s = 0; s < spatial_dimension; ++s) { // *v_s_val += *vel_val * *vel_val; // vel_val++; // } // v_s_val++; // } // Material ** mat_val = &(materials.at(0)); // const Mesh:: ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); // Mesh::ConnectivityTypeList::const_iterator it; // for(it = type_list.begin(); it != type_list.end(); ++it) { // if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue; // UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); // UInt nb_element = fem->getMesh().getNbElement(*it); // Vector * v_square_el = new Vector(nb_element * nb_quadrature_points, 1, "v_square per element"); // fem->interpolateOnQuadraturePoints(*v_square, *v_square_el, 1, *it); // Real * v_square_el_val = v_square_el->values; // UInt * elem_mat_val = element_material[*it]->values; // for (UInt el = 0; el < nb_element; ++el) { // Real rho = mat_val[elem_mat_val[el]]->getRho(); // for (UInt q = 0; q < nb_quadrature_points; ++q) { // *v_square_el_val *= rho; // v_square_el_val++; // } // } // ekin += fem->integrate(*v_square_el, *it); // delete v_square_el; // } // delete v_square; /// reduction sum over all processors allReduce(&ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl; element_material[*it]->printself(stream, indent + 3); stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl; } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/solid_mechanics_model.hh b/model/solid_mechanics/solid_mechanics_model.hh index 60f864f90..222807a86 100644 --- a/model/solid_mechanics/solid_mechanics_model.hh +++ b/model/solid_mechanics/solid_mechanics_model.hh @@ -1,386 +1,398 @@ /** * @file solid_mechanics_model.hh * @author Nicolas Richart * @date[creation] Thu Jul 22 11:51:06 2010 * @date[last modification] Thu Oct 14 14:00:06 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "material.hh" #include "material_parser.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ namespace akantu { // class Material; class IntegrationScheme2ndOrder; class Contact; class Solver; class SparseMatrix; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; SolidMechanicsModel(UInt spatial_dimension, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// allocate all vectors void initVectors(); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model void initModel(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme void updateResidual(bool need_initialize = true); /// compute the acceleration from the residual void updateAcceleration(); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the implicit solver void initImplicitSolver(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// solve Ku = f // void solve(Vector & solution); void solve(); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance); /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type); /// compute force vector from a function(x,y,z) that describe stresses void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// read a custom material with a keyword and class as template template UInt readCustomMaterial(const std::string & filename, const std::string & keyword); private: /// read properties part of a material file and create the material template Material * readMaterialProperties(std::ifstream & infile, MaterialID mat_id, UInt ¤t_line); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); + /// assemble the mass matrix + void assembleMass(); + + private: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); - /// assemble the lumped mass matrix for local and ghost elements - void assembleMassLumpedRowSum(GhostType ghost_type, const ElementType type); + void assembleMass(GhostType ghost_type); + + /// fill a vector of rho + void computeRho(Vector & rho, + ElementType type, + GhostType ghost_type); - /// assemble the lumped mass matrix for local and ghost elements - void assembleMassLumpedDiagonalScaling(GhostType ghost_type, const ElementType type); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const; inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Vector &); /** * @brief get the SolidMechanicsModel::current_position vector \warn only * consistent after a call to SolidMechanicsModel::updateCurrentPosition */ AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); /** * @brief get the SolidMechanicsModel::increment vector \warn only consistent * if SolidMechanicsModel::setIncrementFlagOn has been called before */ AKANTU_GET_MACRO(Increment, *increment, const Vector &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, const Vector &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Vector &); /** * @brief get the SolidMechanicsModel::acceleration vector, updated by * SolidMechanicsModel::updateAcceleration */ AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Vector &); /// get the SolidMechanicsModel::residual vector, computed by SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, const Vector &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Vector &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get the SolidMechanicsModel::element_material vector corresponding to a given akantu::ElementType AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, Vector &); /// get a particular material inline Material & getMaterial(UInt mat_index); /// give the number of materials inline UInt getNbMaterials() { return materials.size(); }; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the equation number Vector AKANTU_GET_MACRO(EquationNumber, *equation_number, const Vector &); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); inline FEM & getFEMBoundary(std::string name = ""); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Vector * displacement; /// lumped mass array Vector * mass; /// velocities array Vector * velocity; /// accelerations array Vector * acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// position of a dof in the K matrix Vector * equation_number; /// local equation_number to global unordered_map::type local_eq_num_to_global; + /// mass matrix + SparseMatrix * mass_matrix; + + /// velocity damping matrix + SparseMatrix * velocity_damping_matrix; + /// stiffness matrix SparseMatrix * stiffness_matrix; /// materials of all elements ByElementTypeUInt element_material; /// materials of all ghost elements ByElementTypeUInt ghost_element_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Vector * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// the spatial dimension UInt spatial_dimension; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/model/solid_mechanics/solid_mechanics_model_mass.cc b/model/solid_mechanics/solid_mechanics_model_mass.cc index 761a6c9b9..e478e5c26 100644 --- a/model/solid_mechanics/solid_mechanics_model_mass.cc +++ b/model/solid_mechanics/solid_mechanics_model_mass.cc @@ -1,95 +1,141 @@ /** * @file solid_mechanics_model_mass.cc * @author Nicolas Richart * @date Mon Oct 4 15:21:23 2010 * * @brief function handling mass computation * * @section LICENSE * * Copyright (©) 2010-2011 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 "solid_mechanics_model.hh" #include "material.hh" //#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); memset(mass->values, 0, nb_nodes*sizeof(Real)); assembleMassLumped(_not_ghost); assembleMassLumped(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); - Material ** mat_val = &(materials.at(0)); + FEM & fem = getFEM(); - UInt nb_element; - UInt * elem_mat_val; + Vector rho_1(0,1); + + const Mesh::ConnectivityTypeList & type_list = fem.getMesh().getConnectivityTypeList(ghost_type); + Mesh::ConnectivityTypeList::const_iterator it; + for(it = type_list.begin(); it != type_list.end(); ++it) { + if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; + ElementType type = *it; - ByElementTypeReal rho_1; + computeRho(rho_1, type, ghost_type); - const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(ghost_type); + fem.assembleFieldLumped(rho_1, *mass, type, ghost_type); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::assembleMass(GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + MyFEMType & fem = getFEMClass(); + + Vector rho_1(0,1); + UInt nb_element; + + const Mesh::ConnectivityTypeList & type_list = fem.getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { + if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; ElementType type = *it; + + const Vector * shapes; if(ghost_type == _not_ghost) { - nb_element = getFEM().getMesh().getNbElement(type); - elem_mat_val = element_material[type]->values; + nb_element = fem.getMesh().getNbElement(type); + shapes = &(fem.getShapes(type)); } else { - nb_element = getFEM().getMesh().getNbGhostElement(type); - elem_mat_val = ghost_element_material[type]->values; + nb_element = fem.getMesh().getNbGhostElement(type); + shapes = &(fem.getGhostShapes(type)); } - UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type); - rho_1[type] = new Vector(nb_element * nb_quadrature_points, 1, "rho_x_1"); - Real * rho_1_val = rho_1[type]->values; + computeRho(rho_1, type, ghost_type); + } - /// compute @f$ rho @f$ for each nodes of each element - for (UInt el = 0; el < nb_element; ++el) { - Real rho = mat_val[elem_mat_val[el]]->getRho(); /// here rho is constant in an element - for (UInt n = 0; n < nb_quadrature_points; ++n) { - *rho_1_val++ = rho; - } - } + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::computeRho(Vector & rho, + ElementType type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + Material ** mat_val = &(materials.at(0)); + UInt * elem_mat_val; + + UInt nb_element; + + FEM & fem = getFEM(); + + if(ghost_type == _not_ghost) { + nb_element = fem.getMesh().getNbElement(type); + elem_mat_val = element_material[type]->values; + } else { + nb_element = fem.getMesh().getNbGhostElement(type); + elem_mat_val = ghost_element_material[type]->values; } - getFEM().assembleFieldLumped(rho_1, *mass, ghost_type); + UInt nb_quadrature_points = fem.getNbQuadraturePoints(type); - for(it = type_list.begin(); it != type_list.end(); ++it) { - delete rho_1[*it]; + rho.resize(nb_element * nb_quadrature_points); + Real * rho_1_val = rho.values; + + /// compute @f$ rho @f$ for each nodes of each element + for (UInt el = 0; el < nb_element; ++el) { + Real rho = mat_val[elem_mat_val[el]]->getRho(); /// here rho is constant in an element + for (UInt n = 0; n < nb_quadrature_points; ++n) { + *rho_1_val++ = rho; + } } AKANTU_DEBUG_OUT(); } + + __END_AKANTU__ diff --git a/test/test_vector/test_matrix.cc b/test/test_vector/test_matrix.cc index 2736230a3..e594aa217 100644 --- a/test/test_vector/test_matrix.cc +++ b/test/test_vector/test_matrix.cc @@ -1,124 +1,125 @@ /** * @file test_matrix.cc * @author Nicolas Richart * @date Fri Feb 18 22:29:00 2011 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 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 #include /* -------------------------------------------------------------------------- */ #include "aka_types.hh" +#include "aka_types_expression.hh" #include "aka_vector.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { #define n 2 UInt nbm = 10000; Real time; Vector A(nbm, n*n); Vector B(nbm, n*n); Vector C1(nbm, n*n); Vector C2(nbm, n*n); Vector C3(nbm, n*n); Vector C4(nbm, n*n); for (UInt i = 0; i < n*n; ++i) { A.values[i] = drand48(); B.values[i] = drand48(); } for (UInt i = 1; i > nbm; ++i) { memcpy(A.values + i * n * n, A.values, n*n*sizeof(Real)); memcpy(B.values + i * n * n, B.values, n*n*sizeof(Real)); } struct timeval begin, end; /* ------------------------------------------------------------------------ */ gettimeofday(&begin, NULL); Math::matrix_matrix(n,n,n,A,B,C1); gettimeofday(&end, NULL); // Vector::iterator > mitC = C1.begin >(); // std::cout << *mitC << std::endl; //time = (end.tv_sec * 1e3 + end.tv_usec * 1e-3) - (begin.tv_sec * 1e3 + begin.tv_usec * 1e-3); time = (end.tv_sec * 1e6 + end.tv_usec) - (begin.tv_sec * 1e6 + begin.tv_usec); std::cout << "matrix_matrix : " << std::fixed << time/nbm << "us" << std::endl; /* ------------------------------------------------------------------------ */ - Vector::iterator itA = A.begin(n,n); - Vector::iterator itB = B.begin(n,n); - Vector::iterator itC = C2.begin(n,n); + Vector::iterator itA = A.begin(n,n); + Vector::iterator itB = B.begin(n,n); + Vector::iterator itC = C2.begin(n,n); gettimeofday(&begin, NULL); for (UInt i = 0; i < nbm; ++i) { *itC = *itA * *itB; ++itA; ++itB;++itC; } gettimeofday(&end, NULL); // itC = C2.begin(n,n); // std::cout << *itC << std::endl; time = (end.tv_sec * 1e6 + end.tv_usec) - (begin.tv_sec * 1e6 + begin.tv_usec); std::cout << "it Ma * it Ma : " << std::fixed << time/nbm << "us" << std::endl; /* ------------------------------------------------------------------------ */ - Vector::iterator > titA = A.begin >(); - Vector::iterator > titB = B.begin >(); - Vector::iterator > titC = C3.begin >(); + Vector::iterator > titA = A. begin >(); + Vector::iterator > titB = B. begin >(); + Vector::iterator > titC = C3.begin >(); gettimeofday(&begin, NULL); for (UInt i = 0; i < nbm; ++i) { *titC = *titA * *titB; ++titA; ++titB;++titC; } gettimeofday(&end, NULL); // titC = C3.begin >(); // std::cout << *titC << std::endl; time = (end.tv_sec * 1e6 + end.tv_usec) - (begin.tv_sec * 1e6 + begin.tv_usec); std::cout << time/nbm << "us" << std::endl; /* ------------------------------------------------------------------------ */ - Vector::iterator muitA = A.begin(n,n); - Vector::iterator muitB = B.begin(n,n); - Vector::iterator muitC = C4.begin(n,n); + Vector::iterator muitA = A.begin(n,n); + Vector::iterator muitB = B.begin(n,n); + Vector::iterator muitC = C4.begin(n,n); gettimeofday(&begin, NULL); for (UInt i = 0; i < nbm; ++i) { (*muitC).mul(*muitA, *muitB); ++muitA; ++muitB;++muitC; } gettimeofday(&end, NULL); // muitC = C4.begin(n,n); // std::cout << *muitC << std::endl; time = (end.tv_sec * 1e6 + end.tv_usec) - (begin.tv_sec * 1e6 + begin.tv_usec); std::cout << time/nbm << "us" << std::endl; return 0; } diff --git a/test/test_vector/test_vector.cc b/test/test_vector/test_vector.cc index 3c046535e..55b164110 100644 --- a/test/test_vector/test_vector.cc +++ b/test/test_vector/test_vector.cc @@ -1,102 +1,102 @@ /** * @file test_vector.cc * @author Nicolas Richart * @date Tue Jun 29 17:32:23 2010 * * @brief test of the vector class * * @section LICENSE * * Copyright (©) 2010-2011 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 /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "aka_vector.hh" /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { int def_value[3]; def_value[0] = 10; def_value[1] = 20; def_value[2] = 30; akantu::Vector int_vect(10, 3, def_value, "test1"); for (unsigned int i = 5; i < int_vect.getSize(); ++i) { for (unsigned int j = 0; j < int_vect.getNbComponent(); ++j) { int_vect.values[i*int_vect.getNbComponent() + j] = def_value[j]*10; } } std::cerr << int_vect; int new_elem[3]; new_elem[0] = 1; new_elem[1] = 2; new_elem[2] = 3; int_vect.push_back(new_elem); int_vect.push_back(200); int_vect.erase(0); std::cerr << int_vect; akantu::Vector int_vect0(0, 3, def_value, "test2"); std::cerr << int_vect0; int_vect0.push_back(new_elem); std::cerr << int_vect0; akantu::Vector mat_vect(10, 4, 1.); memset(mat_vect.values, 0, 10*4*sizeof(double)); - typedef akantu::RealTMatrix<2,2> ItType; + typedef akantu::types::RealTMatrix<2,2> ItType; typedef akantu::Vector RealVector; RealVector::iterator it; it = mat_vect.begin(); RealVector::iterator end = mat_vect.end(); for (; it != end; ++it) { std::cout << *it << std::endl; } - akantu::Matrix m1(2, 3, 1.); - akantu::Matrix m2(3, 5, 2.); - akantu::Matrix m3; + akantu::types::Matrix m1(2, 3, 1.); + akantu::types::Matrix m2(3, 5, 2.); + akantu::types::Matrix m3; m3 = m1 * m2; std::cout << m1 << m2 << m3; - RealVector::iterator itm; + RealVector::iterator itm; itm = mat_vect.begin(2, 2); - RealVector::iterator endm = mat_vect.end (2, 2); + RealVector::iterator endm = mat_vect.end (2, 2); for (; itm != endm; ++itm) { std::cout << *itm << std::endl; } return EXIT_SUCCESS; }