diff --git a/src/common/aka_math.cc b/src/common/aka_math.cc index e2614adcd..3129d10ca 100644 --- a/src/common/aka_math.cc +++ b/src/common/aka_math.cc @@ -1,234 +1,251 @@ /** * @file aka_math.cc * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author David Simon Kammer * @author Nicolas Richart * @author Leonardo Snozzi * @author Peter Spijker * @author Marco Vocialta * * @date creation: Wed Aug 04 2010 * @date last modification: Fri May 15 2015 * * @brief Implementation of the math toolbox * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * 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" #include "aka_array.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -void Math::matrix_vector(UInt m, UInt n, - const Array & A, - const Array & x, - Array & y, - Real alpha) { +void Math::matrix_vector(UInt m, UInt n, const Array & A, + const Array & x, Array & 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"); + "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."); + "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."); + 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.storage(); Real * x_val = x.storage(); Real * y_val = y.storage(); for (UInt el = 0; el < nb_element; ++el) { 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 Array & A, - const Array & B, - Array & C, - Real alpha) { +void Math::matrix_matrix(UInt m, UInt n, UInt k, const Array & A, + const Array & B, Array & 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"); + "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."); + "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( + 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."); + 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.storage(); Real * B_val = B.storage(); Real * C_val = C.storage(); for (UInt el = 0; el < nb_element; ++el) { 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 Array & A, - const Array & B, - Array & C, - Real alpha) { +void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Array & A, + const Array & B, Array & 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"); + "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."); + "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( + 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."); + 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.storage(); Real * B_val = B.storage(); Real * C_val = C.storage(); for (UInt el = 0; el < nb_element; ++el) { 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::compute_tangents(const Array & normals, Array & tangents) { +void Math::compute_tangents(const Array & normals, + Array & tangents) { AKANTU_DEBUG_IN(); - UInt spatial_dimension = normals.getNbComponent(); + UInt spatial_dimension = normals.getNbComponent(); UInt tangent_components = spatial_dimension * (spatial_dimension - 1); - AKANTU_DEBUG_ASSERT(tangent_components == tangents.getNbComponent(), - "Cannot compute the tangents, the storage array for tangents" - << " does not have the good amount of components."); + AKANTU_DEBUG_ASSERT( + tangent_components == tangents.getNbComponent(), + "Cannot compute the tangents, the storage array for tangents" + << " does not have the good amount of components."); UInt nb_normals = normals.getSize(); tangents.resize(nb_normals); - Real * normal_it = normals .storage(); + Real * normal_it = normals.storage(); Real * tangent_it = tangents.storage(); /// compute first tangent for (UInt q = 0; q < nb_normals; ++q) { /// if normal is orthogonal to xy plane, arbitrarly define tangent - if ( Math::are_float_equal(Math::norm2(normal_it), 0) ) + if (Math::are_float_equal(Math::norm2(normal_it), 0)) tangent_it[0] = 1; else Math::normal2(normal_it, tangent_it); - normal_it += spatial_dimension; + normal_it += spatial_dimension; tangent_it += tangent_components; } /// compute second tangent (3D case) if (spatial_dimension == 3) { - normal_it = normals .storage(); + normal_it = normals.storage(); tangent_it = tangents.storage(); for (UInt q = 0; q < nb_normals; ++q) { Math::normal3(normal_it, tangent_it, tangent_it + spatial_dimension); - normal_it += spatial_dimension; + normal_it += spatial_dimension; tangent_it += tangent_components; } } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +Real Math::reduce(Array & array) { + UInt nb_values = array.getSize(); + if (nb_values == 0) + return 0.; + + UInt nb_values_to_sum = nb_values >> 1; + + std::sort(array.begin(), array.end()); + + // as long as the half is not empty + while (nb_values_to_sum) { + UInt remaining = (nb_values - 2 * nb_values_to_sum); + if (remaining) + array(nb_values - 2) += array(nb_values - 1); + // sum to consecutive values and store the sum in the first half + for (UInt i = 0; i < nb_values_to_sum; ++i) { + array(i) = array(2 * i) + array(2 * i + 1); + } + + nb_values = nb_values_to_sum; + nb_values_to_sum >>= 1; + } + + return array(0); +} __END_AKANTU__ diff --git a/src/common/aka_math.hh b/src/common/aka_math.hh index 7ece3baac..c8fd7f403 100644 --- a/src/common/aka_math.hh +++ b/src/common/aka_math.hh @@ -1,302 +1,306 @@ /** * @file aka_math.hh * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Leonardo Snozzi * @author Peter Spijker * @author Marco Vocialta * * @date creation: Wed Aug 04 2010 * @date last modification: Fri May 15 2015 * * @brief mathematical operations * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * 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" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class Array; class Math { /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Matrix algebra */ /* ------------------------------------------------------------------------ */ /// @f$ y = A*x @f$ static void matrix_vector(UInt m, UInt n, const Array & A, const Array & x, Array & y, Real alpha = 1.); /// @f$ y = A*x @f$ static inline void matrix_vector(UInt m, UInt n, Real * A, Real * x, Real * y, Real alpha = 1.); /// @f$ y = A^t*x @f$ static inline void matrixt_vector(UInt m, UInt n, Real * A, Real * x, Real * y, Real alpha = 1.); /// @f$ C = A*B @f$ static void matrix_matrix(UInt m, UInt n, UInt k, const Array & A, const Array & B, Array & C, Real alpha = 1.); /// @f$ C = A*B^t @f$ static void matrix_matrixt(UInt m, UInt n, UInt k, const Array & A, const Array & B, Array & C, Real alpha = 1.); /// @f$ C = A*B @f$ static inline void matrix_matrix(UInt m, UInt n, UInt k, Real * A, Real * B, Real * C, Real alpha = 1.); /// @f$ C = A^t*B @f$ static inline void matrixt_matrix(UInt m, UInt n, UInt k, Real * A, Real * B, Real * C, Real alpha = 1.); /// @f$ C = A*B^t @f$ static inline void matrix_matrixt(UInt m, UInt n, UInt k, Real * A, Real * B, Real * C, Real alpha = 1.); /// @f$ C = A^t*B^t @f$ static inline void matrixt_matrixt(UInt m, UInt n, UInt k, Real * A, Real * B, Real * C, Real alpha = 1.); template static inline void matMul(UInt m, UInt n, UInt k, Real alpha, Real * A, Real * B, Real beta, Real * C); template static inline void matVectMul(UInt m, UInt n, Real alpha, Real * A, Real * x, Real beta, Real * y); static inline void matrix33_eigenvalues(Real * A, Real * Adiag); static inline void matrix22_eigenvalues(Real * A, Real * Adiag); template static inline void eigenvalues(Real * A, Real * d); /// solve @f$ A x = \Lambda x @f$ and return d and V such as @f$ A V[i:] = d[i] V[i:]@f$ template static void matrixEig(UInt n, T * A, T * d, T * V = NULL); /// determinent of a 2x2 matrix static inline Real det2(const Real * mat); /// determinent of a 3x3 matrix static inline Real det3(const Real * mat); /// determinent of a nxn matrix template static inline Real det(const Real * mat); /// determinent of a nxn matrix template static inline T det(UInt n, const T * mat); /// inverse a nxn matrix template static inline void inv(const Real * mat, Real * inv); /// inverse a nxn matrix template static inline void inv(UInt n, const T * mat, T * inv); /// 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); /// solve A x = b using a LU factorization template static inline void solve(UInt n, const T * A, T * x, const T * b); /// return the double dot product between 2 tensors in 2d static inline Real matrixDoubleDot22(Real * A, Real * B); /// return the double dot product between 2 tensors in 3d static inline Real matrixDoubleDot33(Real * A, Real * B); /// extension of the double dot product to two 2nd order tensor in dimension n static inline Real matrixDoubleDot(UInt n, Real * A, Real * B); /* ------------------------------------------------------------------------ */ /* Array algebra */ /* ------------------------------------------------------------------------ */ /// vector cross product static inline void vectorProduct3(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 norm of a vector static inline Real norm(UInt n, 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); /// return the dot product between 2 vectors static inline Real vectorDot(Real * v1, Real * v2, UInt n); /* ------------------------------------------------------------------------ */ /* Geometry */ /* ------------------------------------------------------------------------ */ /// 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); /// compute the tangents to an array of normal vectors static void compute_tangents(const Array & normals, Array & tangents); /// 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); /// test if two vectors are equal within a given tolerance static inline bool are_vector_equal(UInt n, Real * x, Real * y); #ifdef isnan # error "You probably included which is incompatible with aka_math please use\ or add a \"#undef isnan\" before akantu includes" #endif /// test if a real is a NaN static inline bool isnan(Real x); /// test if the line x and y intersects each other static inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max); /// test if a is in the range [x_min, x_max] static inline bool is_in_range(Real a, Real x_min, Real x_max); static inline Real getTolerance() { return tolerance; }; static inline void setTolerance(Real tol) { tolerance = tol; }; template static inline T pow(T x); + /// reduce all the values of an array, the summation is done in place and the + /// array is modified + static Real reduce(Array & array); + class NewtonRaphson { public: NewtonRaphson(Real tolerance, Real max_iteration) : tolerance(tolerance), max_iteration(max_iteration) { } template Real solve(const Functor & funct, Real x_0); private: Real tolerance; Real max_iteration; }; struct NewtonRaphsonFunctor { NewtonRaphsonFunctor(const std::string & name) : name(name) {} virtual Real f(Real x) const = 0; virtual Real f_prime(Real x) const = 0; std::string name; }; private: /// tolerance for functions that need one static Real tolerance; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "aka_math_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_AKA_MATH_H__ */ diff --git a/src/fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc b/src/fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc index f8a8bb084..990d8bc3f 100644 --- a/src/fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc +++ b/src/fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc @@ -1,189 +1,188 @@ /** * @file element_class_quadrangle_8_inline_impl.cc * * @author Nicolas Richart * * @date creation: Wed May 18 2011 * @date last modification: Tue Apr 07 2015 * * @brief Specialization of the ElementClass for the _quadrangle_8 * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * 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 * * @verbatim - \eta - ^ - | + \eta + ^ + | (-1,1) (0,1) (1,1) x-------x-------x | | | | | | | | | (-1,0)| | |(1,0) ----x---------------X-----> \xi | | | | | | | | | | | | x-------x-------x (-1,-1) (0,-1) (1,-1) - | + | @endverbatim * * @subsection shapes Shape functions * @f[ * \begin{array}{lll} * N1 = (1 - \xi) (1 - \eta)(- 1 - \xi - \eta) / 4 * & \frac{\partial N1}{\partial \xi} = (1 - \eta)(2 \xi + \eta) / 4 * & \frac{\partial N1}{\partial \eta} = (1 - \xi)(\xi + 2 \eta) / 4 \\ * N2 = (1 + \xi) (1 - \eta)(- 1 + \xi - \eta) / 4 \\ * & \frac{\partial N2}{\partial \xi} = (1 - \eta)(2 \xi - \eta) / 4 * & \frac{\partial N2}{\partial \eta} = - (1 + \xi)(\xi - 2 \eta) / 4 \\ * N3 = (1 + \xi) (1 + \eta)(- 1 + \xi + \eta) / 4 \\ * & \frac{\partial N3}{\partial \xi} = (1 + \eta)(2 \xi + \eta) / 4 * & \frac{\partial N3}{\partial \eta} = (1 + \xi)(\xi + 2 \eta) / 4 \\ * N4 = (1 - \xi) (1 + \eta)(- 1 - \xi + \eta) / 4 * & \frac{\partial N4}{\partial \xi} = (1 + \eta)(2 \xi - \eta) / 4 * & \frac{\partial N4}{\partial \eta} = - (1 - \xi)(\xi - 2 \eta) / 4 \\ * N5 = (1 - \xi^2) (1 - \eta) / 2 * & \frac{\partial N1}{\partial \xi} = - \xi (1 - \eta) * & \frac{\partial N1}{\partial \eta} = - (1 - \xi^2) / 2 \\ * N6 = (1 + \xi) (1 - \eta^2) / 2 \\ * & \frac{\partial N2}{\partial \xi} = (1 - \eta^2) / 2 * & \frac{\partial N2}{\partial \eta} = - \eta (1 + \xi) \\ * N7 = (1 - \xi^2) (1 + \eta) / 2 \\ * & \frac{\partial N3}{\partial \xi} = - \xi (1 + \eta) * & \frac{\partial N3}{\partial \eta} = (1 - \xi^2) / 2 \\ * N8 = (1 - \xi) (1 - \eta^2) / 2 * & \frac{\partial N4}{\partial \xi} = - (1 - \eta^2) / 2 * & \frac{\partial N4}{\partial \eta} = - \eta (1 - \xi) \\ * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& 0 \qquad \eta_{q0} = 0 * @f} */ /* -------------------------------------------------------------------------- */ -AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_quadrangle_8, _gt_quadrangle_8, _itp_serendip_quadrangle_8, _ek_regular, 2, - _git_segment, 2); +AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_quadrangle_8, _gt_quadrangle_8, + _itp_serendip_quadrangle_8, _ek_regular, 2, + _git_segment, 3); AKANTU_DEFINE_SHAPE(_gt_quadrangle_8, _gst_square); /* -------------------------------------------------------------------------- */ template <> template -inline void -InterpolationElement<_itp_serendip_quadrangle_8>::computeShapes(const vector_type & c, - vector_type & N) { +inline void InterpolationElement<_itp_serendip_quadrangle_8>::computeShapes( + const vector_type & c, vector_type & N) { /// Natural coordinates - const Real xi = c(0); + const Real xi = c(0); const Real eta = c(1); - N(0) = .25 * (1 - xi) * (1 - eta) * (- 1 - xi - eta); - N(1) = .25 * (1 + xi) * (1 - eta) * (- 1 + xi - eta); - N(2) = .25 * (1 + xi) * (1 + eta) * (- 1 + xi + eta); - N(3) = .25 * (1 - xi) * (1 + eta) * (- 1 - xi + eta); - N(4) = .5 * (1 - xi * xi) * (1 - eta ); - N(5) = .5 * (1 + xi ) * (1 - eta * eta); - N(6) = .5 * (1 - xi * xi) * (1 + eta ); - N(7) = .5 * (1 - xi ) * (1 - eta * eta); + N(0) = .25 * (1 - xi) * (1 - eta) * (-1 - xi - eta); + N(1) = .25 * (1 + xi) * (1 - eta) * (-1 + xi - eta); + N(2) = .25 * (1 + xi) * (1 + eta) * (-1 + xi + eta); + N(3) = .25 * (1 - xi) * (1 + eta) * (-1 - xi + eta); + N(4) = .5 * (1 - xi * xi) * (1 - eta); + N(5) = .5 * (1 + xi) * (1 - eta * eta); + N(6) = .5 * (1 - xi * xi) * (1 + eta); + N(7) = .5 * (1 - xi) * (1 - eta * eta); } /* -------------------------------------------------------------------------- */ template <> template -inline void -InterpolationElement<_itp_serendip_quadrangle_8>::computeDNDS(const vector_type & c, - matrix_type & dnds) { +inline void InterpolationElement<_itp_serendip_quadrangle_8>::computeDNDS( + const vector_type & c, matrix_type & dnds) { - const Real xi = c(0); + const Real xi = c(0); const Real eta = c(1); /// dN/dxi dnds(0, 0) = .25 * (1 - eta) * (2 * xi + eta); dnds(0, 1) = .25 * (1 - eta) * (2 * xi - eta); dnds(0, 2) = .25 * (1 + eta) * (2 * xi + eta); dnds(0, 3) = .25 * (1 + eta) * (2 * xi - eta); - dnds(0, 4) = - xi * (1 - eta); - dnds(0, 5) = .5 * (1 - eta * eta); - dnds(0, 6) = - xi * (1 + eta); - dnds(0, 7) = - .5 * (1 - eta * eta); + dnds(0, 4) = -xi * (1 - eta); + dnds(0, 5) = .5 * (1 - eta * eta); + dnds(0, 6) = -xi * (1 + eta); + dnds(0, 7) = -.5 * (1 - eta * eta); /// dN/deta dnds(1, 0) = .25 * (1 - xi) * (2 * eta + xi); dnds(1, 1) = .25 * (1 + xi) * (2 * eta - xi); dnds(1, 2) = .25 * (1 + xi) * (2 * eta + xi); dnds(1, 3) = .25 * (1 - xi) * (2 * eta - xi); - dnds(1, 4) = - .5 * (1 - xi * xi); - dnds(1, 5) = - eta * (1 + xi); - dnds(1, 6) = .5 * (1 - xi * xi); - dnds(1, 7) = - eta * (1 - xi); + dnds(1, 4) = -.5 * (1 - xi * xi); + dnds(1, 5) = -eta * (1 + xi); + dnds(1, 6) = .5 * (1 - xi * xi); + dnds(1, 7) = -eta * (1 - xi); } /* -------------------------------------------------------------------------- */ -template<> +template <> inline Real GeometricalElement<_gt_quadrangle_8>::getInradius(const Matrix & coord) { Real a, b, h; Vector u0 = coord(0); Vector u1 = coord(1); Vector u2 = coord(2); Vector u3 = coord(3); Vector u4 = coord(4); Vector u5 = coord(5); Vector u6 = coord(6); Vector u7 = coord(7); a = u0.distance(u4); b = u4.distance(u1); h = std::min(a, b); a = u1.distance(u5); b = u5.distance(u2); h = std::min(h, std::min(a, b)); a = u2.distance(u6); b = u6.distance(u3); h = std::min(h, std::min(a, b)); a = u3.distance(u7); b = u7.distance(u0); h = std::min(h, std::min(a, b)); return h; } /* -------------------------------------------------------------------------- */ template <> inline void -InterpolationElement<_itp_serendip_quadrangle_8>::computeSpecialJacobian(const Matrix & J, - Real & jac){ +InterpolationElement<_itp_serendip_quadrangle_8>::computeSpecialJacobian( + const Matrix & J, Real & jac) { Vector vprod(J.cols()); Matrix Jt(J.transpose(), true); vprod.crossProduct(Jt(0), Jt(1)); jac = vprod.norm(); } diff --git a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.cc b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.cc index 90e9056de..0ea08bf68 100644 --- a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.cc +++ b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.cc @@ -1,107 +1,107 @@ /** * @file element_class_segment_2_inline_impl.cc * * @author Nicolas Richart * * @date creation: Fri Jul 16 2010 * @date last modification: Sun Oct 19 2014 * * @brief Specialization of the element_class class for the type _segment_2 * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * 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 * * @verbatim q --x--------|--------x---> x -1 0 1 @endverbatim * * @subsection shapes Shape functions * @f{eqnarray*}{ * w_1(x) &=& 1/2(1 - x) \\ * w_2(x) &=& 1/2(1 + x) * @f} * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * x_{q} &=& 0 * @f} */ /* -------------------------------------------------------------------------- */ -AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_segment_2, _gt_segment_2, _itp_lagrange_segment_2, _ek_regular, 1, - _git_segment, 1); +AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_segment_2, _gt_segment_2, + _itp_lagrange_segment_2, _ek_regular, 1, + _git_segment, 1); AKANTU_DEFINE_SHAPE(_gt_segment_2, _gst_square); /* -------------------------------------------------------------------------- */ template <> template -inline void -InterpolationElement<_itp_lagrange_segment_2>::computeShapes(const vector_type & natural_coords, - vector_type & N) { +inline void InterpolationElement<_itp_lagrange_segment_2>::computeShapes( + const vector_type & natural_coords, vector_type & N) { /// natural coordinate Real c = natural_coords(0); /// shape functions - N(0) = 0.5*(1-c); - N(1) = 0.5*(1+c); + N(0) = 0.5 * (1 - c); + N(1) = 0.5 * (1 + c); } /* -------------------------------------------------------------------------- */ template <> template -inline void -InterpolationElement<_itp_lagrange_segment_2>::computeDNDS(__attribute__ ((unused)) const vector_type & natural_coords, - matrix_type & dnds){ - +inline void InterpolationElement<_itp_lagrange_segment_2>::computeDNDS( + __attribute__((unused)) const vector_type & natural_coords, + matrix_type & dnds) { + /// dN1/de - dnds(0, 0) = - .5; + dnds(0, 0) = -.5; /// dN2/de - dnds(0, 1) = .5; + dnds(0, 1) = .5; } /* -------------------------------------------------------------------------- */ template <> inline void -InterpolationElement<_itp_lagrange_segment_2>::computeSpecialJacobian(const Matrix & dxds, - Real & jac) { +InterpolationElement<_itp_lagrange_segment_2>::computeSpecialJacobian( + const Matrix & dxds, Real & jac) { jac = dxds.norm(); } /* -------------------------------------------------------------------------- */ -template<> +template <> inline Real GeometricalElement<_gt_segment_2>::getInradius(const Matrix & coord) { return std::abs(coord(0, 0) - coord(0, 1)); } -// /* -------------------------------------------------------------------------- */ -// template<> inline bool ElementClass<_segment_2>::contains(const Vector & natural_coords) { +// /* -------------------------------------------------------------------------- +// */ +// template<> inline bool ElementClass<_segment_2>::contains(const Vector +// & natural_coords) { // if (natural_coords(0) < -1.) return false; // if (natural_coords(0) > 1.) return false; // return true; // } /* -------------------------------------------------------------------------- */ - - diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh index 16093b382..d16a486e2 100644 --- a/src/fe_engine/fe_engine_template_tmpl.hh +++ b/src/fe_engine/fe_engine_template_tmpl.hh @@ -1,1669 +1,1668 @@ /** * @file fe_engine_template_tmpl.hh * * @author Guillaume Anciaux * @author Mauro Corrado * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Nov 19 2015 * * @brief Template implementation of FEEngineTemplate * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template