diff --git a/common/aka_math.hh b/common/aka_math.hh index c37a7aa20..517d3268d 100644 --- a/common/aka_math.hh +++ b/common/aka_math.hh @@ -1,163 +1,163 @@ /** * @file aka_math.h * @author Nicolas Richart * @date Wed Jul 28 11:51:56 2010 * * @brief mathematical operations * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #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); /// @f$ y = A*x @f$ static inline void matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y); /// @f$ C = A*B @f$ static void matrix_matrix(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C); /// @f$ C = A*B @f$ static inline void matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C); /// @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); /// @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); /// @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); /// 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 * coord); + 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 * coord); + 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); }; /* -------------------------------------------------------------------------- */ /* 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 a3dde0bb1..ae01931d5 100644 --- a/common/aka_math_inline_impl.cc +++ b/common/aka_math_inline_impl.cc @@ -1,368 +1,373 @@ /** * @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 * * \ * */ #ifdef AKANTU_USE_CBLAS # ifndef AKANTU_USE_CBLAS_MKL # include # else // AKANTU_USE_CBLAS_MKL # include # endif //AKANTU_USE_CBLAS_MKL #endif //AKANTU_USE_CBLAS /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y) { #ifdef AKANTU_USE_CBLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1, 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]; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1, 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]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, 1, 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]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, 1, 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]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, m, n, k, 1, 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]; } } } #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 * coord) { - Real xx[9],vol; - const Real *x1 = coord,*x2 = coord+3,*x3 = coord+6,*x4 = coord+9; - - xx[0] = x2[0];xx[1] = x2[1];xx[2] = x2[2]; - xx[3] = x3[0];xx[4] = x3[1];xx[5] = x3[2]; - xx[6] = x4[0];xx[7] = x4[1];xx[8] = x4[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] = x1[0];xx[1] = x1[1];xx[2] = x1[2]; - xx[3] = x3[0];xx[4] = x3[1];xx[5] = x3[2]; - xx[6] = x4[0];xx[7] = x4[1];xx[8] = x4[2]; + 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] = x1[0];xx[1] = x1[1];xx[2] = x1[2]; - xx[3] = x2[0];xx[4] = x2[1];xx[5] = x2[2]; - xx[6] = x4[0];xx[7] = x4[1];xx[8] = x4[2]; + 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] = x1[0];xx[1] = x1[1];xx[2] = x1[2]; - xx[3] = x2[0];xx[4] = x2[1];xx[5] = x2[2]; - xx[6] = x3[0];xx[7] = x3[1];xx[8] = x3[2]; + 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 * coord) { +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(coord , coord+3); - l13 = distance_3d(coord , coord+6); - l14 = distance_3d(coord , coord+9); - l23 = distance_3d(coord+3, coord+6); - l24 = distance_3d(coord+3, coord+9); - l34 = distance_3d(coord+6, coord+9); + 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(coord); + 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[3] = y[2]-x[2]; } diff --git a/fem/element_class.cc b/fem/element_class.cc index 36c6f631a..a155a29e1 100644 --- a/fem/element_class.cc +++ b/fem/element_class.cc @@ -1,124 +1,127 @@ /** * @file element_class.cc * @author Nicolas Richart * @date Tue Jul 20 10:12:44 2010 * * @brief Common part of element_classes * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template UInt ElementClass::nb_nodes_per_element = 0; template ElementType ElementClass::p1_element_type = _not_defined; template UInt ElementClass::nb_quadrature_points = 0; template UInt ElementClass::spatial_dimension = 0; template ElementType ElementClass::facet_type = _not_defined; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_point>::nb_nodes_per_element = 1; template<> ElementType ElementClass<_point>::p1_element_type = _point; template<> ElementType ElementClass<_point>::facet_type = _not_defined; template<> UInt ElementClass<_point>::spatial_dimension = 0; template<> UInt ElementClass<_point>::nb_facets = 0; template<> UInt * ElementClass<_point>::facet_connectivity[] = {}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_segment_2>::nb_nodes_per_element = 2; template<> ElementType ElementClass<_segment_2>::p1_element_type = _segment_2; template<> UInt ElementClass<_segment_2>::nb_quadrature_points = 1; template<> Real ElementClass<_segment_2>::quad[] = {0}; template<> UInt ElementClass<_segment_2>::spatial_dimension = 1; template<> UInt ElementClass<_segment_2>::nb_facets = 2; template<> ElementType ElementClass<_segment_2>::facet_type = _point; template<> UInt ElementClass<_segment_2>::vec_facet_connectivity[]= {0, 1}; template<> UInt * ElementClass<_segment_2>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[1]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_segment_3>::nb_nodes_per_element = 3; template<> ElementType ElementClass<_segment_3>::p1_element_type = _segment_2; template<> UInt ElementClass<_segment_3>::nb_quadrature_points = 2; template<> Real ElementClass<_segment_3>::quad[] = {-1./sqrt(3.),1./sqrt(3.)}; template<> UInt ElementClass<_segment_3>::spatial_dimension = 1; template<> UInt ElementClass<_segment_3>::nb_facets = 2; template<> ElementType ElementClass<_segment_3>::facet_type = _point; template<> UInt ElementClass<_segment_3>::vec_facet_connectivity[]= {0, 1}; template<> UInt * ElementClass<_segment_3>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[1]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_3>::nb_nodes_per_element = 3; template<> ElementType ElementClass<_triangle_3>::p1_element_type = _triangle_3; template<> UInt ElementClass<_triangle_3>::nb_quadrature_points = 1; template<> Real ElementClass<_triangle_3>::quad[] = {1./3., 1./3.}; template<> UInt ElementClass<_triangle_3>::spatial_dimension = 2; template<> UInt ElementClass<_triangle_3>::nb_facets = 3; template<> ElementType ElementClass<_triangle_3>::facet_type = _segment_2; template<> UInt ElementClass<_triangle_3>::vec_facet_connectivity[]= {0, 1, 1, 2, 2, 0}; template<> UInt * ElementClass<_triangle_3>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[2], &vec_facet_connectivity[4]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_6>::nb_nodes_per_element = 6; template<> ElementType ElementClass<_triangle_6>::p1_element_type = _triangle_3; template<> UInt ElementClass<_triangle_6>::nb_quadrature_points = 3; template<> Real ElementClass<_triangle_6>::quad[] = {1./6., 1./6., 2./3., 1./6., 1./6., 2./3.}; template<> UInt ElementClass<_triangle_6>::spatial_dimension = 2; template<> UInt ElementClass<_triangle_6>::nb_facets = 3; template<> ElementType ElementClass<_triangle_6>::facet_type = _segment_3; template<> UInt ElementClass<_triangle_6>::vec_facet_connectivity[]= {0, 1, 3, 1, 2, 4, 2, 0, 5}; template<> UInt * ElementClass<_triangle_6>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_4>::nb_nodes_per_element = 4; template<> ElementType ElementClass<_tetrahedron_4>::p1_element_type = _tetrahedron_4; template<> UInt ElementClass<_tetrahedron_4>::nb_quadrature_points = 1; -template<> Real ElementClass<_tetrahedron_4>::quad[] = {1./4., 1./4.,1./4.}; +template<> Real ElementClass<_tetrahedron_4>::quad[] = {1./4., 1./4., 1./4.}; template<> UInt ElementClass<_tetrahedron_4>::spatial_dimension = 3; template<> UInt ElementClass<_tetrahedron_4>::nb_facets = 4; template<> ElementType ElementClass<_tetrahedron_4>::facet_type = _triangle_3; template<> UInt ElementClass<_tetrahedron_4>::vec_facet_connectivity[]= {0, 2, 1, 1, 2, 3, 2, 0, 3, 0, 1, 3}; template<> UInt * ElementClass<_tetrahedron_4>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6], &vec_facet_connectivity[9]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_10>::nb_nodes_per_element = 10; template<> ElementType ElementClass<_tetrahedron_10>::p1_element_type = _tetrahedron_4; template<> UInt ElementClass<_tetrahedron_10>::nb_quadrature_points = 4; -template<> Real ElementClass<_tetrahedron_10>::quad[] = {}; +template<> Real ElementClass<_tetrahedron_10>::quad[] = {0.1381966011250, 0.1381966011250, 0.1381966011250, // a = (5-sqrt(5))/20 + 0.5854101966250, 0.1381966011250, 0.1381966011250, // b = (5+3*sqrt(5))/20 + 0.1381966011250, 0.5854101966250, 0.1381966011250, + 0.1381966011250, 0.1381966011250, 0.5854101966250}; template<> UInt ElementClass<_tetrahedron_10>::spatial_dimension = 3; template<> UInt ElementClass<_tetrahedron_10>::nb_facets = 4; template<> ElementType ElementClass<_tetrahedron_10>::facet_type = _triangle_6; template<> UInt ElementClass<_tetrahedron_10>::vec_facet_connectivity[]= {0, 2, 1, 6, 5, 4, - 1, 2, 3, 5, 8, 9, - 2, 0, 3, 8, 7, 6, - 0, 1, 3, 4, 9, 7}; + 1, 2, 3, 5, 9, 8, + 2, 0, 3, 6, 7, 9, + 0, 1, 3, 4, 8, 7}; template<> UInt * ElementClass<_tetrahedron_10>::facet_connectivity[] = {&vec_facet_connectivity[0], - &vec_facet_connectivity[3], &vec_facet_connectivity[6], - &vec_facet_connectivity[9]}; + &vec_facet_connectivity[12], + &vec_facet_connectivity[18]}; /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/fem/element_class_inline_impl.cc b/fem/element_class_inline_impl.cc index 8dd7a18d4..52a9d6982 100644 --- a/fem/element_class_inline_impl.cc +++ b/fem/element_class_inline_impl.cc @@ -1,222 +1,222 @@ /** * @file element_class_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 15 10:28:28 2010 * * @brief Implementation of the inline functions of the class element_class * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getNbQuadraturePoint() { return nb_quadrature_points; } /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getShapeSize() { return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getShapeDerivativesSize() { return nb_nodes_per_element * spatial_dimension; } /* -------------------------------------------------------------------------- */ template void ElementClass::preComputeStandards(const Real * coord, const UInt dimension, Real * shape, Real * dshape, Real * jacobians) { // ask for computation of shapes computeShapes(quad, nb_quadrature_points, shape); // compute dnds Real dnds[nb_nodes_per_element * spatial_dimension * nb_quadrature_points]; computeDNDS(quad, nb_quadrature_points, dnds); // compute dxds Real dxds[dimension * spatial_dimension * nb_quadrature_points]; computeDXDS(dnds, nb_quadrature_points, coord, dimension, dxds); // jacobian computeJacobian(dxds, nb_quadrature_points, dimension, jacobians); // if dimension == spatial_dimension compute shape derivatives if (dimension == spatial_dimension) { computeShapeDerivatives(dxds, dnds, nb_quadrature_points, dimension, dshape); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes) { Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapes(cpoint, shapes); shapes += nb_nodes_per_element; cpoint += spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(const Real * natural_coords, const UInt nb_points, Real * dnds) { Real * cpoint = const_cast(natural_coords); Real * cdnds = dnds; for (UInt p = 0; p < nb_points; ++p) { computeDNDS(cpoint, cdnds); cpoint += spatial_dimension; cdnds += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const UInt nb_points, const Real * node_coords, const UInt dimension, Real * dxds) { Real * cdnds = const_cast(dnds); Real * cdxds = dxds; for (UInt p = 0; p < nb_points; ++p) { computeDXDS(cdnds, node_coords, dimension, cdxds); cdnds += nb_nodes_per_element * spatial_dimension; cdxds += spatial_dimension * dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const Real * node_coords, const UInt dimension, Real * dxds) { /// @f$ J = dxds = dnds * x @f$ Math::matrix_matrix(spatial_dimension, dimension, nb_nodes_per_element, dnds, node_coords, dxds); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Real * dxds, const UInt nb_points, const UInt dimension, Real * jac) { Real * cdxds = const_cast(dxds); Real * cjac = jac; for (UInt p = 0; p < nb_points; ++p) { computeJacobian(cdxds, dimension, *cjac); AKANTU_DEBUG_ASSERT((cjac[0] > 0), "Negative jacobian computed, possible problem in the element node order."); cdxds += spatial_dimension * dimension; cjac++; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, const UInt nb_points, const UInt dimension, Real * shape_deriv) { AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"gradient in space " << dimension << " cannot be evaluated for element of dimension " << spatial_dimension); Real * cdxds = const_cast(dxds); Real * cdnds = const_cast(dnds); for (UInt p = 0; p < nb_points; ++p) { computeShapeDerivatives(cdxds, cdnds, shape_deriv); cdnds += spatial_dimension * nb_nodes_per_element; cdxds += spatial_dimension * spatial_dimension; shape_deriv += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, Real * shape_deriv) { - /// @f$ dxds = J^{-1} @f$ + /// @f$ dxds = J^{-1} @f$ Real inv_dxds[spatial_dimension * spatial_dimension]; if (spatial_dimension == 1) inv_dxds[0] = 1./dxds[0]; if (spatial_dimension == 2) Math::inv2(dxds, inv_dxds); if (spatial_dimension == 3) Math::inv3(dxds, inv_dxds); Math::matrixt_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, dnds, inv_dxds, shape_deriv); } /* -------------------------------------------------------------------------- */ template inline Real ElementClass::getInradius(__attribute__ ((unused)) const Real * coord) { AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); return 0; } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeNormalsOnQuadPoint(const Real * coord, const UInt dimension, Real * normals) { AKANTU_DEBUG_ASSERT((dimension - 1) == spatial_dimension, "cannot extract a normal because of dimension mismatch " << dimension << " " << spatial_dimension); Real * cpoint = const_cast(quad); Real * cnormals = normals; Real dnds[spatial_dimension*nb_nodes_per_element]; Real dxds[spatial_dimension*dimension]; for (UInt p = 0; p < nb_quadrature_points; ++p) { computeDNDS(cpoint,dnds); computeDXDS(dnds,coord,dimension,dxds); if (dimension == 2) { Math::normal2(dxds,cnormals); } if (dimension == 3){ Math::normal3(dxds,dxds+dimension,cnormals); } cpoint += spatial_dimension; cnormals += dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * shapes){ AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * dnds){ AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(__attribute__ ((unused)) const Real * dxds, __attribute__ ((unused)) const UInt dimension, __attribute__ ((unused)) Real & jac){ AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } #include "element_classes/element_class_segment_2.cc" #include "element_classes/element_class_segment_3.cc" #include "element_classes/element_class_triangle_3.cc" #include "element_classes/element_class_triangle_6.cc" #include "element_classes/element_class_tetrahedron_4.cc" -// #include "element_classes/element_class_tetrahedron_10.cc" +#include "element_classes/element_class_tetrahedron_10.cc" diff --git a/fem/element_classes/element_class_tetrahedron_10.cc b/fem/element_classes/element_class_tetrahedron_10.cc index 6a0c74199..22fcaf457 100644 --- a/fem/element_classes/element_class_tetrahedron_10.cc +++ b/fem/element_classes/element_class_tetrahedron_10.cc @@ -1,118 +1,255 @@ /** - * @file element_class_inline_impl.cc - * @author Nicolas Richart - * @date Thu Jul 15 10:28:28 2010 + * @file element_class_tetrahedron_10.cc + * @author Peter Spijker + * @date Thu Dec 1 10:28:28 2010 * * @brief Specialization of the element_class class for the type _tetrahedron_10 * * @section LICENSE * * \ * * @section DESCRIPTION * * @verbatim \zeta ^ | (0,0,1) x |` . | ` . | ` . | ` . (0,0.5,0.5) | ` x. | q4 o ` . \eta | ` . -, (0,0,0.5) x ` x (0.5,0,0.5) - | ` x-(0,1,0) | q3 o` - ' - | 0,0.5,0) - ` ' + | (0,0.5,0) - ` ' | x- ` x (0.5,0.5,0) | q1 o - o q2` ' | - ` ' | - ` ' x---------------x--------------` x-----> \xi (0,0,0) (0.5,0,0) (1,0,0) @endverbatim * * @subsection coords Nodes coordinates * * @f[ * \begin{array}{lll} * \xi_{0} = 0 & \eta_{0} = 0 & \zeta_{0} = 0 \\ * \xi_{1} = 1 & \eta_{1} = 0 & \zeta_{1} = 0 \\ * \xi_{2} = 0 & \eta_{2} = 1 & \zeta_{2} = 0 \\ * \xi_{3} = 0 & \eta_{3} = 0 & \zeta_{3} = 1 \\ * \xi_{4} = 1/2 & \eta_{4} = 0 & \zeta_{4} = 0 \\ * \xi_{5} = 1/2 & \eta_{5} = 1/2 & \zeta_{5} = 0 \\ * \xi_{6} = 0 & \eta_{6} = 1/2 & \zeta_{6} = 0 \\ * \xi_{7} = 0 & \eta_{7} = 0 & \zeta_{7} = 1/2 \\ * \xi_{8} = 1/2 & \eta_{8} = 0 & \zeta_{8} = 1/2 \\ * \xi_{9} = 0 & \eta_{9} = 1/2 & \zeta_{9} = 1/2 * \end{array} * @f] * * @subsection shapes Shape functions * @f[ * \begin{array}{llll} * N1 = (1 - \xi - \eta - \zeta) (1 - 2 \xi - 2 \eta - 2 \zeta) - * & \frac{\partial N1}{\partial \xi} = - * & \frac{\partial N1}{\partial \eta} = - * & \frac{\partial N1}{\partial \zeta} = \\ + * & \frac{\partial N1}{\partial \xi} = 4 \xi + 4 \eta + 4 \zeta - 3 + * & \frac{\partial N1}{\partial \eta} = 4 \xi + 4 \eta + 4 \zeta - 3 + * & \frac{\partial N1}{\partial \zeta} = 4 \xi + 4 \eta + 4 \zeta - 3 \\ * N2 = \xi (2 \xi - 1) - * & \frac{\partial N2}{\partial \xi} = + * & \frac{\partial N2}{\partial \xi} = 4 \xi - 1 * & \frac{\partial N2}{\partial \eta} = 0 * & \frac{\partial N2}{\partial \zeta} = 0 \\ * N3 = \eta (2 \eta - 1) * & \frac{\partial N3}{\partial \xi} = 0 - * & \frac{\partial N3}{\partial \eta} = + * & \frac{\partial N3}{\partial \eta} = 4 \eta - 1 * & \frac{\partial N3}{\partial \zeta} = 0 \\ * N4 = \zeta (2 \zeta - 1) * & \frac{\partial N4}{\partial \xi} = 0 * & \frac{\partial N4}{\partial \eta} = 0 - * & \frac{\partial N4}{\partial \zeta} = \\ + * & \frac{\partial N4}{\partial \zeta} = 4 \zeta - 1 \\ * N5 = 4 \xi (1 - \xi - \eta - \zeta) - * & \frac{\partial N5}{\partial \xi} = + * & \frac{\partial N5}{\partial \xi} = 4 - 8 \xi - 4 \eta - 4 \zeta * & \frac{\partial N5}{\partial \eta} = -4 \xi * & \frac{\partial N5}{\partial \zeta} = -4 \xi \\ * N6 = 4 \xi \eta * & \frac{\partial N6}{\partial \xi} = 4 \eta * & \frac{\partial N6}{\partial \eta} = 4 \xi * & \frac{\partial N6}{\partial \zeta} = 0 \\ * N7 = 4 \eta (1 - \xi - \eta - \zeta) - * & \frac{\partial N7}{\partial \xi} = -4 \xi - * & \frac{\partial N7}{\partial \eta} = - * & \frac{\partial N7}{\partial \zeta} = -4 \xi \\ + * & \frac{\partial N7}{\partial \xi} = -4 \eta + * & \frac{\partial N7}{\partial \eta} = 4 - 4 \xi - 8 \eta - 4 \zeta + * & \frac{\partial N7}{\partial \zeta} = -4 \eta \\ * N8 = 4 \zeta (1 - \xi - \eta - \zeta) * & \frac{\partial N8}{\partial \xi} = -4 \zeta * & \frac{\partial N8}{\partial \eta} = -4 \zeta - * & \frac{\partial N8}{\partial \zeta} = \\ + * & \frac{\partial N8}{\partial \zeta} = 4 - 4 \xi - 4 \eta - 8 \zeta \\ * N9 = 4 \zeta \xi * & \frac{\partial N9}{\partial \xi} = 4 \zeta * & \frac{\partial N9}{\partial \eta} = 0 * & \frac{\partial N9}{\partial \zeta} = 4 \xi \\ * N10 = 4 \eta \zeta * & \frac{\partial N10}{\partial \xi} = 0 * & \frac{\partial N10}{\partial \eta} = 4 \zeta * & \frac{\partial N10}{\partial \zeta} = 4 \eta \\ * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f[ * a = \frac{5 - \sqrt{5}}{20}\\ * b = \frac{5 + 3 \sqrt{5}}{20} * \begin{array}{lll} * \xi_{q_0} = a & \eta_{q_0} = a & \zeta_{q_0} = a \\ * \xi_{q_1} = b & \eta_{q_1} = a & \zeta_{q_1} = a \\ * \xi_{q_2} = a & \eta_{q_2} = b & \zeta_{q_2} = a \\ * \xi_{q_3} = a & \eta_{q_3} = a & \zeta_{q_3} = b * \end{array} * @f] */ /* -------------------------------------------------------------------------- */ +template<> UInt ElementClass<_tetrahedron_10>::nb_nodes_per_element; +template<> UInt ElementClass<_tetrahedron_10>::nb_quadrature_points; +template<> UInt ElementClass<_tetrahedron_10>::spatial_dimension; +/* -------------------------------------------------------------------------- */ +template <> inline void ElementClass<_tetrahedron_10>::computeShapes(const Real * natural_coords, + Real * shapes){ + /// Natural coordinates + Real xi = natural_coords[0]; + Real eta = natural_coords[1]; + Real zeta = natural_coords[2]; + Real sum = xi + eta + zeta; + Real c0 = 1 - sum; + Real c1 = 1 - 2*sum; + Real c2 = 2*xi - 1; + Real c3 = 2*eta - 1; + Real c4 = 2*zeta - 1; + + /// Shape functions + shapes[0] = c0 * c1; + shapes[1] = xi * c2; + shapes[2] = eta * c3; + shapes[3] = zeta * c4; + shapes[4] = 4 * xi * c0; + shapes[5] = 4 * xi * eta; + shapes[6] = 4 * eta * c0; + shapes[7] = 4 * zeta * c0; + shapes[8] = 4 * xi * zeta; + shapes[9] = 4 * eta * zeta; +} /* -------------------------------------------------------------------------- */ +template <> inline void ElementClass<_tetrahedron_10>::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, + Real * dnds) { + + /** + * @f[ + * dnds = \left( + * \begin{array}{cccccccccc} + * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial \xi} + * & \frac{\partial N3}{\partial \xi} & \frac{\partial N4}{\partial \xi} + * & \frac{\partial N5}{\partial \xi} & \frac{\partial N6}{\partial \xi} + * & \frac{\partial N7}{\partial \xi} & \frac{\partial N8}{\partial \xi} + * & \frac{\partial N9}{\partial \xi} & \frac{\partial N10}{\partial \xi} \\ + * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} + * & \frac{\partial N3}{\partial \eta} & \frac{\partial N4}{\partial \eta} + * & \frac{\partial N5}{\partial \eta} & \frac{\partial N6}{\partial \eta} + * & \frac{\partial N7}{\partial \eta} & \frac{\partial N8}{\partial \eta} + * & \frac{\partial N9}{\partial \eta} & \frac{\partial N10}{\partial \eta} \\ + * \frac{\partial N1}{\partial \zeta} & \frac{\partial N2}{\partial \zeta} + * & \frac{\partial N3}{\partial \zeta} & \frac{\partial N4}{\partial \zeta} + * & \frac{\partial N5}{\partial \zeta} & \frac{\partial N6}{\partial \zeta} + * & \frac{\partial N7}{\partial \zeta} & \frac{\partial N8}{\partial \zeta} + * & \frac{\partial N9}{\partial \zeta} & \frac{\partial N10}{\partial \zeta} + * \end{array} + * \right) + * @f] + */ + + /// Natural coordinates + Real xi = natural_coords[0]; + Real eta = natural_coords[1]; + Real zeta = natural_coords[2]; + Real sum = xi + eta + zeta; + + /// dN/dxi + dnds[0] = 4 * sum - 3; + dnds[1] = 4 * xi - 1; + dnds[2] = 0; + dnds[3] = 0; + dnds[4] = 4 * (1 - sum - xi); + dnds[5] = 4 * eta; + dnds[6] = -4 * eta; + dnds[7] = -4 * zeta; + dnds[8] = 4 * zeta; + dnds[9] = 0; + + /// dN/deta + dnds[10] = 4 * sum - 3; + dnds[11] = 0; + dnds[12] = 4 * eta - 1; + dnds[13] = 0; + dnds[14] = -4 * xi; + dnds[15] = 4 * xi; + dnds[16] = 4 * (1 - sum - eta); + dnds[17] = -4 * zeta; + dnds[18] = 0; + dnds[19] = 4 * zeta; + + /// dN/dzeta + dnds[20] = 4 * sum - 3; + dnds[21] = 0; + dnds[22] = 0; + dnds[23] = 4 * zeta - 1; + dnds[24] = -4 * xi; + dnds[25] = 0; + dnds[26] = -4 * eta; + dnds[27] = 4 * (1 - sum - zeta); + dnds[28] = 4 * xi; + dnds[29] = 4 * eta; + +} + +/* -------------------------------------------------------------------------- */ +template <> inline void ElementClass<_tetrahedron_10>::computeJacobian(const Real * dxds, + const UInt dimension, + Real & jac) { + + if (dimension == spatial_dimension){ + Real weight = 1./24.; + Real det_dxds = Math::det3(dxds); + jac = det_dxds * weight; + } + else { + AKANTU_DEBUG_ERROR("to be implemented"); + } +} + +/* -------------------------------------------------------------------------- */ +template<> inline Real ElementClass<_tetrahedron_10>::getInradius(const Real * coord) { + + // Only take the four corner tetrahedra + UInt tetrahedra[4][4] = { + {0, 4, 6, 7}, + {4, 1, 5, 8}, + {6, 5, 2, 9}, + {7, 8, 9, 3} + }; + + Real inradius = std::numeric_limits::max(); + for (UInt t = 0; t < 4; t++) { + Real ir = Math::tetrahedron_inradius(coord + tetrahedra[t][0] * spatial_dimension, + coord + tetrahedra[t][1] * spatial_dimension, + coord + tetrahedra[t][2] * spatial_dimension, + coord + tetrahedra[t][3] * spatial_dimension); + inradius = ir < inradius ? ir : inradius; + } + + return inradius; +} diff --git a/fem/element_classes/element_class_tetrahedron_4.cc b/fem/element_classes/element_class_tetrahedron_4.cc index 4839bf081..355741e7c 100644 --- a/fem/element_classes/element_class_tetrahedron_4.cc +++ b/fem/element_classes/element_class_tetrahedron_4.cc @@ -1,118 +1,118 @@ /** * @file element_class_tetrahedron_4.cc * @author Guillaume ANCIAUX * @date Mon Aug 16 18:09:53 2010 * * @brief Specialization of the element_class class for the type _tetrahedron_4 * * @section LICENSE * * \ * * @section DESCRIPTION * * @verbatim \eta ^ | x (0,0,1,0) |` | ` ° \xi | ` ° - | ` x (0,0,0,1) | q.` - ' | -` ' | - ` ' | - ` ' x------------------x-----> \zeta (1,0,0,0) (0,1,0,0) @endverbatim * * @subsection shapes Shape functions * @f{eqnarray*}{ * N1 &=& 1 - \xi - \eta - \zeta \\ * N2 &=& \xi \\ * N3 &=& \eta \\ * N4 &=& \zeta * @f} * * @subsection quad_points Position of quadrature points * @f[ * \xi_{q0} = 1/4 \qquad \eta_{q0} = 1/4 \qquad \zeta_{q0} = 1/4 * @f] */ /* -------------------------------------------------------------------------- */ // /// shape functions // shape[0] = 1./4.; /// N1(q_0) // shape[1] = 1./4.; /// N2(q_0) // shape[2] = 1./4.; /// N3(q_0) // shape[3] = 1./4.; /// N4(q_0) /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_4>::nb_nodes_per_element; template<> UInt ElementClass<_tetrahedron_4>::nb_quadrature_points; template<> UInt ElementClass<_tetrahedron_4>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_4>::computeShapes(const Real * natural_coords, Real * shapes){ Real c0 = natural_coords[1]; /// @f$ c0 = \eta @f$ Real c1 = natural_coords[2]; /// @f$ c1 = \zeta @f$ Real c2 = 1 - natural_coords[0] - natural_coords[1] - natural_coords[2];/// @f$ c2 = 1 - \xi - \eta - \zeta @f$ Real c3 = natural_coords[0]; /// @f$ c2 = \xi @f$ shapes[0] = c0; shapes[1] = c1; shapes[2] = c2; shapes[3] = c3; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_4>::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} & \frac{\partial N4}{\partial \xi} \\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} & \frac{\partial N4}{\partial \eta} \\ * \frac{\partial N1}{\partial \zeta} & \frac{\partial N2}{\partial \zeta} * & \frac{\partial N3}{\partial \zeta} & \frac{\partial N4}{\partial \zeta} * \end{array} * \right) * @f] */ dnds[0] = -1.; dnds[1] = 1.; dnds[2] = 0.; dnds[3] = 0.; dnds[4] = -1.; dnds[5] = 0.; dnds[6] = 1.; dnds[7] = 0.; dnds[8] = -1.; dnds[9] = 0.; dnds[10] = 0.; dnds[11] = 1.; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_4>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac) { if (dimension == spatial_dimension){ Real weight = 1./6.; Real det_dxds = Math::det3(dxds); jac = det_dxds * weight; } else { AKANTU_DEBUG_ERROR("to be implemented"); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_tetrahedron_4>::getInradius(const Real * coord) { - return Math::tetrahedron_inradius(coord); + return Math::tetrahedron_inradius(coord, coord+3, coord+6, coord+9); } diff --git a/fem/element_classes/element_class_triangle_6.cc b/fem/element_classes/element_class_triangle_6.cc index 6500ff09f..ee6cc6fc9 100644 --- a/fem/element_classes/element_class_triangle_6.cc +++ b/fem/element_classes/element_class_triangle_6.cc @@ -1,193 +1,193 @@ /** * @file element_class_triangle_6.cc * @author Nicolas Richart * @date Thu Jul 15 10:28:28 2010 * * @brief Specialization of the element_class class for the type _triangle_6 * * @section LICENSE * * \ * * @section DESCRIPTION * * @verbatim \eta ^ | x 2 | ` | ` | . ` | q2 ` 5 x x 4 | ` | ` | .q0 q1. ` | ` x---------x---------x-----> \xi 0 3 1 @endverbatim * * @subsection coords Nodes coordinates * * @f[ * \begin{array}{ll} * \xi_{0} = 0 & \eta_{0} = 0 \\ * \xi_{1} = 1 & \eta_{1} = 0 \\ * \xi_{2} = 0 & \eta_{2} = 1 \\ * \xi_{3} = 1/2 & \eta_{3} = 0 \\ * \xi_{4} = 1/2 & \eta_{4} = 1/2 \\ * \xi_{5} = 0 & \eta_{5} = 1/2 * \end{array} * @f] * * @subsection shapes Shape functions * @f[ * \begin{array}{lll} * N1 = -(1 - \xi - \eta) (1 - 2 (1 - \xi - \eta)) * & \frac{\partial N1}{\partial \xi} = 1 - 4(1 - \xi - \eta) * & \frac{\partial N1}{\partial \eta} = 1 - 4(1 - \xi - \eta) \\ * N2 = - \xi (1 - 2 \xi) * & \frac{\partial N1}{\partial \xi} = - 1 + 4 \xi * & \frac{\partial N1}{\partial \eta} = 0 \\ * N3 = - \eta (1 - 2 \eta) * & \frac{\partial N1}{\partial \xi} = 0 * & \frac{\partial N1}{\partial \eta} = - 1 + 4 \eta \\ * N4 = 4 \xi (1 - \xi - \eta) * & \frac{\partial N1}{\partial \xi} = 4 (1 - 2 \xi - \eta) * & \frac{\partial N1}{\partial \eta} = - 4 \eta \\ * N5 = 4 \xi \eta * & \frac{\partial N1}{\partial \xi} = 4 \xi * & \frac{\partial N1}{\partial \eta} = 4 \eta \\ * N6 = 4 \eta (1 - \xi - \eta) * & \frac{\partial N1}{\partial \xi} = - 4 \xi * & \frac{\partial N1}{\partial \eta} = 4 (1 - \xi - 2 \eta) * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& 1/6 \qquad \eta_{q0} = 1/6 \\ * \xi_{q1} &=& 2/3 \qquad \eta_{q1} = 1/6 \\ * \xi_{q2} &=& 1/6 \qquad \eta_{q2} = 2/3 * @f} */ // /// quadrature point position // quad[0] = 1./6.; /// q0_{\xi} // quad[1] = 1./6.; /// q0_{\eta} // quad[2] = 2./3.; /// q1_{\xi} // quad[3] = 1./6.; /// q1_{\eta} // quad[4] = 1./6.; /// q2_{\xi} // quad[5] = 2./3.; /// q2_{\eta} /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_6>::nb_nodes_per_element; template<> UInt ElementClass<_triangle_6>::nb_quadrature_points; template<> UInt ElementClass<_triangle_6>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_6>::computeShapes(const Real * natural_coords, Real * shapes){ /// Natural coordinates Real c0 = 1 - natural_coords[0] - natural_coords[1]; /// @f$ c0 = 1 - \xi - \eta @f$ Real c1 = natural_coords[0]; /// @f$ c1 = \xi @f$ Real c2 = natural_coords[1]; /// @f$ c2 = \eta @f$ shapes[0] = c0 * (2 * c0 - 1.); shapes[1] = c1 * (2 * c1 - 1.); shapes[2] = c2 * (2 * c2 - 1.); shapes[3] = 4 * c0 * c1; shapes[4] = 4 * c1 * c2; shapes[5] = 4 * c2 * c0; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_6>::computeDNDS(const Real * natural_coords, Real * dnds){ /** * @f[ * dnds = \left( * \begin{array}{cccccc} * \frac{\partial N1}{\partial \xi} * & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} * & \frac{\partial N4}{\partial \xi} * & \frac{\partial N5}{\partial \xi} * & \frac{\partial N6}{\partial \xi} \ \ * * \frac{\partial N1}{\partial \eta} * & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} * & \frac{\partial N4}{\partial \eta} * & \frac{\partial N5}{\partial \eta} * & \frac{\partial N6}{\partial \eta} * \end{array} * \right) @f] */ /// Natural coordinates Real c0 = 1 - natural_coords[0] - natural_coords[1]; /// @f$ c0 = 1 - \xi - \eta @f$ Real c1 = natural_coords[0]; /// @f$ c1 = \xi @f$ Real c2 = natural_coords[1]; /// @f$ c2 = \eta @f$ dnds[0] = 1 - 4 * c0; dnds[1] = 4 * c1 - 1.; dnds[2] = 0.; dnds[3] = 4 * (c0 - c1); dnds[4] = 4 * c2; dnds[5] = - 4 * c2; dnds[6] = 1 - 4 * c0; dnds[7] = 0.; dnds[8] = 4 * c2 - 1.; dnds[9] = - 4 * c1; dnds[10] = 4 * c1; dnds[11] = 4 * (c0 - c2); } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_6>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ // if element dimension is the same as the space dimension // then jacobian factor is the determinent of dxds if (dimension == spatial_dimension){ Real weight = 1./6.; jac = Math::det2(dxds)*weight; AKANTU_DEBUG_ASSERT(jac > 0, "Negative jacobian computed, possible problem in the element node order."); } else { AKANTU_DEBUG_ERROR("to implement"); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_triangle_6>::getInradius(const Real * coord) { UInt triangles[4][3] = { {0, 3, 5}, {3, 1, 4}, {3, 4, 5}, {5, 4, 2} }; Real inradius = std::numeric_limits::max(); for (UInt t = 0; t < 4; t++) { - Real ir = Math::triangle_inradius(coord + triangles[t][0] * 2, - coord + triangles[t][1] * 2, - coord + triangles[t][2] * 2); + Real ir = Math::triangle_inradius(coord + triangles[t][0] * spatial_dimension, + coord + triangles[t][1] * spatial_dimension, + coord + triangles[t][2] * spatial_dimension); inradius = ir < inradius ? ir : inradius; } return inradius; } diff --git a/model/solid_mechanics/solid_mechanics_model_mass.cc b/model/solid_mechanics/solid_mechanics_model_mass.cc index 64b68e7bd..1fa1f2105 100644 --- a/model/solid_mechanics/solid_mechanics_model_mass.cc +++ b/model/solid_mechanics/solid_mechanics_model_mass.cc @@ -1,201 +1,206 @@ /** * @file solid_mechanics_model_mass.cc * @author Nicolas Richart * @date Mon Oct 4 15:21:23 2010 * * @brief function handling mass computation * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "material.hh" //#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = fem->getMesh().getNbNodes(); memset(mass->values, 0, nb_nodes*sizeof(Real)); assembleMassLumped(_not_ghost); assembleMassLumped(_ghost); /// for not connected nodes put mass to one in order to avoid /// wrong range in paraview Real * mass_values = mass->values; for (UInt i = 0; i < nb_nodes; ++i) { if (fabs(mass_values[i]) < std::numeric_limits::epsilon() || isnan(mass_values[i])) mass_values[i] = 1; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); 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; switch(*it) { case _triangle_6: + case _tetrahedron_10: assembleMassLumpedDiagonalScaling(ghost_type, *it); break; default: assembleMassLumpedRowSum(ghost_type, *it); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$ */ void SolidMechanicsModel::assembleMassLumpedRowSum(GhostType ghost_type, ElementType type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt shape_size = FEM::getShapeSize(type); UInt nb_element; const Vector * shapes; UInt * elem_mat_val; if(ghost_type == _not_ghost) { nb_element = fem->getMesh().getNbElement(type); shapes = &(fem->getShapes(type)); elem_mat_val = element_material[type]->values; } else { nb_element = fem->getMesh().getNbGhostElement(type); shapes = &(fem->getGhostShapes(type)); elem_mat_val = ghost_element_material[type]->values; } if (nb_element == 0) { AKANTU_DEBUG_OUT(); return; } Vector * rho_phi_i = new Vector(nb_element * nb_quadrature_points, shape_size, "rho_x_shapes"); Real * rho_phi_i_val = rho_phi_i->values; Real * shapes_val = shapes->values; /// compute @f$ rho * \varphi_i @f$ for each nodes of each element for (UInt el = 0; el < nb_element; ++el) { Real rho = mat_val[elem_mat_val[el]]->getRho(); for (UInt n = 0; n < shape_size * nb_quadrature_points; ++n) { *rho_phi_i_val++ = rho * *shapes_val++; } } Vector * int_rho_phi_i = new Vector(0, shape_size, "inte_rho_x_shapes"); fem->integrate(*rho_phi_i, *int_rho_phi_i, shape_size, type, ghost_type); delete rho_phi_i; fem->assembleVector(*int_rho_phi_i, *mass, 1, type, ghost_type); delete int_rho_phi_i; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$ */ void SolidMechanicsModel::assembleMassLumpedDiagonalScaling(GhostType ghost_type, ElementType type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); UInt nb_nodes_per_element_p1 = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type)); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt nb_element; UInt * elem_mat_val; 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(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; } if (nb_element == 0) { AKANTU_DEBUG_OUT(); return; } Vector * rho_1 = new Vector(nb_element * nb_quadrature_points, 1, "rho_x_1"); Real * rho_1_val = rho_1->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; } } /// compute @f$ \int \rho dV = \rho V @f$ for each element Vector * int_rho_1 = new Vector(nb_element * nb_quadrature_points, 1, "inte_rho_x_1"); fem->integrate(*rho_1, *int_rho_1, 1, type, ghost_type); delete rho_1; /// distribute the mass of the element to the nodes Vector * mass_per_node = new Vector(nb_element, nb_nodes_per_element, "mass_per_node"); Real * int_rho_1_val = int_rho_1->values; Real * mass_per_node_val = mass_per_node->values; for (UInt e = 0; e < nb_element; ++e) { Real lmass = *int_rho_1_val * corner_factor; for (UInt n = 0; n < nb_nodes_per_element_p1; ++n) *mass_per_node_val++ = lmass; /// corner points lmass = *int_rho_1_val * mid_factor; for (UInt n = nb_nodes_per_element_p1; n < nb_nodes_per_element; ++n) *mass_per_node_val++ = lmass; /// mid points int_rho_1_val++; } delete int_rho_1; fem->assembleVector(*mass_per_node, *mass, 1, type, ghost_type); delete mass_per_node; AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/test_fem/CMakeLists.txt b/test/test_fem/CMakeLists.txt index 76a152f29..b7abef2a9 100644 --- a/test/test_fem/CMakeLists.txt +++ b/test/test_fem/CMakeLists.txt @@ -1,90 +1,102 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # 1D #=============================================================================== add_mesh(test_fem_line_1_mesh line.geo 1 1 OUTPUT line1.msh) add_mesh(test_fem_line_2_mesh line.geo 1 2 OUTPUT line2.msh) #=============================================================================== # 1st order register_test(test_integrate_segment_2 test_integrate_segment_2.cc) add_dependencies(test_integrate_segment_2 test_fem_line_1_mesh) #=============================================================================== register_test(test_interpolate_segment_2 test_interpolate_segment_2.cc) add_dependencies(test_interpolate_segment_2 test_fem_line_1_mesh) #=============================================================================== register_test(test_gradient_segment_2 test_gradient_segment_2.cc) add_dependencies(test_gradient_segment_2 test_fem_line_1_mesh) #=============================================================================== # 2nd order register_test(test_integrate_segment_3 test_integrate_segment_3.cc) add_dependencies(test_integrate_segment_3 test_fem_line_2_mesh) #=============================================================================== register_test(test_interpolate_segment_3 test_interpolate_segment_3.cc) add_dependencies(test_interpolate_segment_3 test_fem_line_2_mesh) #=============================================================================== register_test(test_gradient_segment_3 test_gradient_segment_3.cc) add_dependencies(test_gradient_segment_3 test_fem_line_2_mesh) #=============================================================================== # 2D #=============================================================================== add_mesh(test_fem_square_1_mesh square.geo 2 1 OUTPUT square1.msh) add_mesh(test_fem_square_2_mesh square.geo 2 2 OUTPUT square2.msh) add_mesh(test_fem_circle_1_mesh circle.geo 2 1 OUTPUT circle1.msh) add_mesh(test_fem_circle_2_mesh circle.geo 2 2 OUTPUT circle2.msh) #=============================================================================== # 1st order register_test(test_integrate_triangle_3 test_integrate_triangle_3.cc) add_dependencies(test_integrate_triangle_3 test_fem_square_1_mesh test_fem_circle_1_mesh) #=============================================================================== register_test(test_interpolate_triangle_3 test_interpolate_triangle_3.cc) add_dependencies(test_interpolate_triangle_3 test_fem_square_1_mesh) #=============================================================================== register_test(test_gradient_triangle_3 test_gradient_triangle_3.cc) add_dependencies(test_gradient_triangle_3 test_fem_square_1_mesh) #=============================================================================== # 2nd order register_test(test_integrate_triangle_6 test_integrate_triangle_6.cc) add_dependencies(test_integrate_triangle_6 test_fem_square_2_mesh test_fem_circle_2_mesh) #=============================================================================== register_test(test_interpolate_triangle_6 test_interpolate_triangle_6.cc) add_dependencies(test_interpolate_triangle_6 test_fem_square_2_mesh) #=============================================================================== register_test(test_gradient_triangle_6 test_gradient_triangle_6.cc) add_dependencies(test_gradient_triangle_6 test_fem_square_2_mesh) #=============================================================================== # 3D #=============================================================================== -add_mesh(test_fem_cube_2_mesh cube.geo 3 1 OUTPUT cube1.msh) +add_mesh(test_fem_cube_1_mesh cube.geo 3 1 OUTPUT cube1.msh) +add_mesh(test_fem_cube_2_mesh cube.geo 3 2 OUTPUT cube2.msh) #=============================================================================== # 1st order register_test(test_integrate_tetrahedron_4 test_integrate_tetrahedron_4.cc) -add_dependencies(test_integrate_tetrahedron_4 test_fem_cube_2_mesh) +add_dependencies(test_integrate_tetrahedron_4 test_fem_cube_1_mesh) #=============================================================================== register_test(test_interpolate_tetrahedron_4 test_interpolate_tetrahedron_4.cc) -add_dependencies(test_interpolate_tetrahedron_4 test_fem_cube_2_mesh) +add_dependencies(test_interpolate_tetrahedron_4 test_fem_cube_1_mesh) #=============================================================================== register_test(test_gradient_tetrahedron_4 test_gradient_tetrahedron_4.cc) -add_dependencies(test_gradient_tetrahedron_4 test_fem_cube_2_mesh) +add_dependencies(test_gradient_tetrahedron_4 test_fem_cube_1_mesh) + +#=============================================================================== +# 2nd order +register_test(test_integrate_tetrahedron_10 test_integrate_tetrahedron_10.cc) +add_dependencies(test_integrate_tetrahedron_10 test_fem_cube_2_mesh) +#=============================================================================== +register_test(test_interpolate_tetrahedron_10 test_interpolate_tetrahedron_10.cc) +add_dependencies(test_interpolate_tetrahedron_10 test_fem_cube_2_mesh) +#=============================================================================== +register_test(test_gradient_tetrahedron_10 test_gradient_tetrahedron_10.cc) +add_dependencies(test_gradient_tetrahedron_10 test_fem_cube_2_mesh) diff --git a/test/test_fem/cube.geo b/test/test_fem/cube.geo index 3c98f9918..1fa3a481f 100644 --- a/test/test_fem/cube.geo +++ b/test/test_fem/cube.geo @@ -1,37 +1,37 @@ -h = 0.5; +h = 0.25; Point(1) = {0.0, 0.0, 0.0, h}; Point(2) = {1.0, 0.0, 0.0, h}; Point(3) = {0.0, 1.0, 0.0, h}; Point(4) = {1.0, 1.0, 0.0, h}; Point(5) = {0.0, 0.0, 1.0, h}; Point(6) = {1.0, 0.0, 1.0, h}; Point(7) = {0.0, 1.0, 1.0, h}; Point(8) = {1.0, 1.0, 1.0, h}; Line(1) = {7, 8}; Line(2) = {8, 4}; Line(3) = {4, 3}; Line(4) = {3, 7}; Line(5) = {1, 5}; Line(6) = {5, 6}; Line(7) = {6, 2}; Line(8) = {2, 1}; Line(9) = {3, 1}; Line(10) = {7, 5}; Line(11) = {8, 6}; Line(12) = {4, 2}; Line Loop(13) = {1, 11, -6, -10}; Plane Surface(14) = {13}; Line Loop(15) = {3, 4, 1, 2}; Plane Surface(16) = {15}; Line Loop(17) = {6, 7, 8, 5}; Plane Surface(18) = {17}; Line Loop(19) = {3, 9, -8, -12}; Plane Surface(20) = {19}; Line Loop(21) = {4, 10, -5, -9}; Plane Surface(22) = {21}; Line Loop(23) = {11, 7, -12, -2}; Plane Surface(24) = {23}; Surface Loop(25) = {24, 14, 16, 20, 22, 18}; Volume(26) = {25}; diff --git a/test/test_fem/test_gradient_tetrahedron_10.cc b/test/test_fem/test_gradient_tetrahedron_10.cc new file mode 100644 index 000000000..325b378ef --- /dev/null +++ b/test/test_fem/test_gradient_tetrahedron_10.cc @@ -0,0 +1,68 @@ +/** + * @file test_gradient_tetrahedron_4.cc + * @author Nicolas Richart + * @date Mon Jul 19 10:55:49 2010 + * + * @brief test of the fem class + * + * @section LICENSE + * + * \ + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "fem.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" + + +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + UInt dim = 3; + ElementType type = _tetrahedron_10; + MeshIOMSH mesh_io; + Mesh my_mesh(dim); + mesh_io.read("cube2.msh", my_mesh); + FEM *fem = new FEM(my_mesh, dim, "my_fem"); + + debug::setDebugLevel(dblDump); + fem->initShapeFunctions(); + + std::cout << *fem << std::endl; + + StaticMemory * st_mem = StaticMemory::getStaticMemory(); + std::cout << *st_mem << std::endl; + + Vector const_val(fem->getMesh().getNbNodes(), 2, "const_val"); + Vector grad_on_quad(0, 2 * dim, "grad_on_quad"); + + for (UInt i = 0; i < const_val.getSize(); ++i) { + const_val.values[i * 2 + 0] = 1.; + const_val.values[i * 2 + 1] = 2.; + } + + fem->gradientOnQuadraturePoints(const_val, grad_on_quad, 2, type); + std::ofstream my_file("out.txt"); + my_file << const_val << std::endl; + my_file << grad_on_quad << std::endl; + + // compute gradient of coordinates + Vector grad_coord_on_quad(0, 9, "grad_coord_on_quad"); + fem->gradientOnQuadraturePoints(my_mesh.getNodes(), grad_coord_on_quad, my_mesh.getSpatialDimension(), type); + my_file << my_mesh.getNodes() << std::endl; + my_file << grad_coord_on_quad << std::endl; + + delete fem; + finalize(); + + return EXIT_SUCCESS; +} diff --git a/test/test_fem/test_integrate_tetrahedron_10.cc b/test/test_fem/test_integrate_tetrahedron_10.cc new file mode 100644 index 000000000..9702067c5 --- /dev/null +++ b/test/test_fem/test_integrate_tetrahedron_10.cc @@ -0,0 +1,72 @@ +/** + * @file test_integrate_tetrahedron_4.cc + * @author Nicolas Richart + * @date Mon Jul 19 10:55:49 2010 + * + * @brief test of the fem class + * + * @section LICENSE + * + * \ + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "fem.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" + + +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + UInt dim = 3; + ElementType type = _tetrahedron_10; + MeshIOMSH mesh_io; + Mesh my_mesh(dim); + mesh_io.read("cube2.msh", my_mesh); + FEM *fem = new FEM(my_mesh, dim, "my_fem"); + + debug::_debug_level = dblDump; + fem->initShapeFunctions(); + + std::cout << *fem << std::endl; + + StaticMemory * st_mem = StaticMemory::getStaticMemory(); + std::cout << *st_mem << std::endl; + + Vector const_val(fem->getMesh().getNbNodes(), 2, "const_val"); + Vector val_on_quad(0, 2, "val_on_quad"); + + for (UInt i = 0; i < const_val.getSize(); ++i) { + const_val.values[i * 2 + 0] = 1.; + const_val.values[i * 2 + 1] = 2.; + } + //interpolate function on quadrature points + fem->interpolateOnQuadraturePoints(const_val, val_on_quad, 2, type); + //integrate function on elements + akantu::Vector int_val_on_elem(0, 2, "int_val_on_elem"); + fem->integrate(val_on_quad, int_val_on_elem, 2, type); + // get global integration value + Real value[2] = {0,0}; + std::ofstream my_file("out.txt"); + my_file << int_val_on_elem << std::endl; + for (UInt i = 0; i < fem->getMesh().getNbElement(type); ++i) { + value[0] += int_val_on_elem.values[2*i]; + value[1] += int_val_on_elem.values[2*i+1]; + } + + my_file << "integral is " << value[0] << " " << value[1] << std::endl; + + delete fem; + finalize(); + + return EXIT_SUCCESS; +} diff --git a/test/test_fem/test_interpolate_tetrahedron_10.cc b/test/test_fem/test_interpolate_tetrahedron_10.cc new file mode 100644 index 000000000..af52c9a10 --- /dev/null +++ b/test/test_fem/test_interpolate_tetrahedron_10.cc @@ -0,0 +1,69 @@ +/** + * @file test_interpolate_tetrahedron_4.cc + * @author Nicolas Richart + * @date Mon Jul 19 10:55:49 2010 + * + * @brief test of the fem class + * + * @section LICENSE + * + * \ + * + */ + +/* -------------------------------------------------------------------------- */ +#include +#include +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "fem.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" + + +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + UInt dim = 3; + ElementType type = _tetrahedron_10; + MeshIOMSH mesh_io; + Mesh my_mesh(dim); + mesh_io.read("cube2.msh", my_mesh); + FEM *fem = new FEM(my_mesh, dim, "my_fem"); + + debug::setDebugLevel(dblDump); + fem->initShapeFunctions(); + + std::cout << *fem << std::endl; + + StaticMemory * st_mem = StaticMemory::getStaticMemory(); + std::cout << *st_mem << std::endl; + + Vector const_val(fem->getMesh().getNbNodes(), dim, "const_val"); + Vector val_on_quad(0, dim, "val_on_quad"); + + for (UInt i = 0; i < const_val.getSize(); ++i) { + for (UInt j = 0; j < dim; ++j) { + const_val.values[i * dim + j] = j; + } + } + + fem->interpolateOnQuadraturePoints(const_val, val_on_quad, dim, type); + std::ofstream my_file("out.txt"); + my_file << const_val << std::endl; + my_file << val_on_quad << std::endl; + + // interpolate coordinates + Vector coord_on_quad(0, my_mesh.getSpatialDimension(), "coord_on_quad"); + fem->interpolateOnQuadraturePoints(my_mesh.getNodes(), coord_on_quad, my_mesh.getSpatialDimension(), type); + my_file << my_mesh.getNodes() << std::endl; + my_file << coord_on_quad << std::endl; + + delete fem; + finalize(); + + return EXIT_SUCCESS; +} diff --git a/test/test_model/test_solid_mechanics_model/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/CMakeLists.txt index 5d032a295..c5c721ce0 100644 --- a/test/test_model/test_solid_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/CMakeLists.txt @@ -1,63 +1,68 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== add_subdirectory("test_materials") #=============================================================================== add_mesh(test_solid_mechanics_model_square_mesh triangle.geo 2 1) add_mesh(test_solid_mechanics_model_circle_mesh circle.geo 2 2) register_test(test_solid_mechanics_model_square test_solid_mechanics_model_square.cc) add_dependencies(test_solid_mechanics_model_square test_solid_mechanics_model_square_mesh test_solid_mechanics_model_circle_mesh) register_test(test_solid_mechanics_model_circle_2 test_solid_mechanics_model_circle_2.cc) #=============================================================================== add_mesh(test_bar_traction_2d_mesh1 bar.geo 2 1 OUTPUT bar1.msh) add_mesh(test_bar_traction_2d_mesh2 bar.geo 2 2 OUTPUT bar2.msh) register_test(test_solid_mechanics_model_bar_traction2d test_solid_mechanics_model_bar_traction2d.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d test_bar_traction_2d_mesh1 test_bar_traction_2d_mesh2) if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) register_test(test_solid_mechanics_model_bar_traction2d_parallel test_solid_mechanics_model_bar_traction2d_parallel.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d_parallel test_bar_traction_2d_mesh2) add_mesh(test_solid_mechanics_model_segment_mesh segment.geo 1 2) register_test(test_solid_mechanics_model_segment_parallel test_solid_mechanics_model_segment_parallel.cc) add_dependencies(test_solid_mechanics_model_segment_parallel test_solid_mechanics_model_segment_mesh) endif() #=============================================================================== -add_mesh(test_cube3d_mesh cube.geo 3 1) +add_mesh(test_cube3d_mesh1 cube.geo 3 1 OUTPUT cube1.msh) +add_mesh(test_cube3d_mesh2 cube.geo 3 2 OUTPUT cube2.msh) register_test(test_solid_mechanics_model_cube3d test_solid_mechanics_model_cube3d.cc) -add_dependencies(test_solid_mechanics_model_cube3d test_cube3d_mesh) +add_dependencies(test_solid_mechanics_model_cube3d test_cube3d_mesh1) + +register_test(test_solid_mechanics_model_cube3d_tetra10 + test_solid_mechanics_model_cube3d_tetra10.cc) +add_dependencies(test_solid_mechanics_model_cube3d_tetra10 test_cube3d_mesh2) #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc index b7bb84dfa..f84412fda 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc @@ -1,115 +1,115 @@ /** * @file test_solid_mechanics_model_cube3d.cc * @author Guillaume ANCIAUX * @date Tue Aug 17 11:31:22 2010 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER int main(int argc, char *argv[]) { akantu::UInt max_steps = 10000; akantu::Real epot, ekin; + akantu::ElementType type = akantu::_tetrahedron_4; + akantu::UInt paratype = TETRA1; akantu::Mesh mesh(3); akantu::MeshIOMSH mesh_io; - mesh_io.read("cube.msh", mesh); + mesh_io.read("cube1.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); /// model initialization model->initVectors(); /// initialize the vectors akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, 3*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); akantu::Real time_step = model->getStableTimeStep(); model->setTimeStep(time_step/10.); model->assembleMassLumped(); std::cout << *model << std::endl; /// boundary conditions akantu::Real eps = 1e-16; for (akantu::UInt i = 0; i < nb_nodes; ++i) { model->getDisplacement().values[3*i] = model->getFEM().getMesh().getNodes().values[3*i] / 100.; if(model->getFEM().getMesh().getNodes().values[3*i] <= eps) { model->getBoundary().values[3*i ] = true; - if(model->getFEM().getMesh().getNodes().values[3*i + 1] <= eps) - model->getBoundary().values[3*i + 1] = true; } if(model->getFEM().getMesh().getNodes().values[3*i + 1] <= eps) { model->getBoundary().values[3*i + 1] = true; } } // model->getDisplacement().values[1] = 0.1; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; // dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, 3, nb_nodes, "coordinates"); - dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(akantu::_tetrahedron_4).values, - TETRA1, model->getFEM().getMesh().getNbElement(akantu::_tetrahedron_4), C_MODE); + dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, + paratype, model->getFEM().getMesh().getNbElement(type), C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, 3, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, 3, "velocity"); dumper.AddNodeDataField(model->getMass().values, 1, "mass"); dumper.AddNodeDataField(model->getResidual().values, 3, "force"); - dumper.AddElemDataField(model->getMaterial(0).getStrain(akantu::_tetrahedron_4).values, 9, "strain"); - dumper.AddElemDataField(model->getMaterial(0).getStress(akantu::_tetrahedron_4).values, 9, "stress"); + dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, 9, "strain"); + dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, 9, "stress"); dumper.SetPrefix("paraview/"); dumper.Init(); #endif //AKANTU_USE_IOHELPER for(akantu::UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); std::cout << s << " " << epot << " " << ekin << " " << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc similarity index 79% copy from test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc copy to test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc index b7bb84dfa..9aca0ce17 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc @@ -1,115 +1,115 @@ /** * @file test_solid_mechanics_model_cube3d.cc * @author Guillaume ANCIAUX * @date Tue Aug 17 11:31:22 2010 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @section LICENSE * * \ * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER int main(int argc, char *argv[]) { akantu::UInt max_steps = 10000; akantu::Real epot, ekin; + akantu::ElementType type = akantu::_tetrahedron_10; + akantu::UInt paratype = TETRA2; akantu::Mesh mesh(3); akantu::MeshIOMSH mesh_io; - mesh_io.read("cube.msh", mesh); + mesh_io.read("cube2.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); /// model initialization model->initVectors(); /// initialize the vectors akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, 3*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); akantu::Real time_step = model->getStableTimeStep(); model->setTimeStep(time_step/10.); model->assembleMassLumped(); std::cout << *model << std::endl; /// boundary conditions - akantu::Real eps = 1e-16; + akantu::Real eps = 1e-2; for (akantu::UInt i = 0; i < nb_nodes; ++i) { - model->getDisplacement().values[3*i] = model->getFEM().getMesh().getNodes().values[3*i] / 100.; + model->getDisplacement().values[3*i+2] = model->getFEM().getMesh().getNodes().values[3*i+2] / 100.; - if(model->getFEM().getMesh().getNodes().values[3*i] <= eps) { - model->getBoundary().values[3*i ] = true; - if(model->getFEM().getMesh().getNodes().values[3*i + 1] <= eps) - model->getBoundary().values[3*i + 1] = true; + if(model->getFEM().getMesh().getNodes().values[3*i + 2] <= eps) { + model->getBoundary().values[3*i + 2] = true; } - if(model->getFEM().getMesh().getNodes().values[3*i + 1] <= eps) { - model->getBoundary().values[3*i + 1] = true; + if(model->getFEM().getMesh().getNodes().values[3*i] <= eps) { + model->getBoundary().values[3*i] = true; } } // model->getDisplacement().values[1] = 0.1; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; // dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, 3, nb_nodes, "coordinates"); - dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(akantu::_tetrahedron_4).values, - TETRA1, model->getFEM().getMesh().getNbElement(akantu::_tetrahedron_4), C_MODE); + dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, + paratype, model->getFEM().getMesh().getNbElement(type), C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, 3, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, 3, "velocity"); dumper.AddNodeDataField(model->getMass().values, 1, "mass"); dumper.AddNodeDataField(model->getResidual().values, 3, "force"); - dumper.AddElemDataField(model->getMaterial(0).getStrain(akantu::_tetrahedron_4).values, 9, "strain"); - dumper.AddElemDataField(model->getMaterial(0).getStress(akantu::_tetrahedron_4).values, 9, "stress"); + dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, 9, "strain"); + dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, 9, "stress"); dumper.SetPrefix("paraview/"); dumper.Init(); #endif //AKANTU_USE_IOHELPER for(akantu::UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); std::cout << s << " " << epot << " " << ekin << " " << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } return EXIT_SUCCESS; }