diff --git a/common/aka_common.hh b/common/aka_common.hh index 97da8f2c4..029ecc123 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,347 +1,350 @@ /** * @file aka_common.hh * @author Nicolas Richart * @date Fri Jun 11 09:48:06 2010 * * @namespace akantu * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ } /* -------------------------------------------------------------------------- */ #include "aka_error.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef double Real; typedef unsigned int UInt; typedef unsigned long long UInt64; typedef signed int Int; typedef std::string ID; static const Real UINT_INIT_VALUE = 0; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = 0; #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; typedef ID VectorID; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef ID MeshID; typedef ID FEMID; typedef ID ModelID; typedef ID MaterialID; typedef ID SparseMatrixID; typedef ID SolverID; typedef ID ShapeID; typedef ID IntegratorID; typedef UInt Surface; /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #include #define AKANTU_BOOST_CASE_MACRO(r,macro,type) \ case type : { macro(type); break;} #define AKANTU_BOOST_ELEMENT_SWITCH(macro); \ do { \ switch(type) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO,macro,AKANTU_ELEMENT_TYPE) \ case _not_defined: \ case _max_element_type: { \ AKANTU_DEBUG_ERROR("Wrong type : " << type); \ break; \ } \ } \ } while(0) #define AKANTU_BOOST_LIST_MACRO(r,macro,type) \ macro(type) #define AKANTU_BOOST_ELEMENT_LIST(macro) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO,macro,AKANTU_ELEMENT_TYPE) /* -------------------------------------------------------------------------- */ /// @boost sequence of element to loop on in global tasks #define AKANTU_ELEMENT_TYPE \ (_segment_2) \ (_segment_3) \ (_triangle_3) \ (_triangle_6) \ (_tetrahedron_4) \ (_tetrahedron_10) \ (_quadrangle_4) \ + (_quadrangle_8) \ (_hexahedron_8) \ (_point) /// @enum ElementType type of element potentially contained in a Mesh enum ElementType { _not_defined = 0, _segment_2 = 1, /// first order segment _segment_3 = 2, /// second order segment _triangle_3 = 3, /// first order triangle _triangle_6 = 4, /// second order triangle _tetrahedron_4 = 5, /// first order tetrahedron _tetrahedron_10 = 6, /// second order tetrahedron @remark not implemented yet _quadrangle_4, /// first order quadrangle + _quadrangle_8, /// second order quadrangle _hexahedron_8, /// first order hexahedron _point, /// point only for some algorithm to be generic like mesh partitioning _max_element_type }; /// @enum MaterialType different materials implemented enum MaterialType { _elastic = 0, _max_material_type }; typedef void (*BoundaryFunction)(double *,double *); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_forces }; /// @enum SparseMatrixType type of sparse matrix used enum SparseMatrixType { _unsymmetric, _symmetric }; /* -------------------------------------------------------------------------- */ /* Contact */ /* -------------------------------------------------------------------------- */ typedef ID ContactID; typedef ID ContactSearchID; typedef ID ContactNeighborStructureID; enum ContactType { _ct_not_defined = 0, _ct_2d_expli = 1, _ct_3d_expli = 2, _ct_rigid = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum GhostSynchronizationTag type of synchronizations enum GhostSynchronizationTag { /// SolidMechanicsModel tags _gst_smm_mass, /// synchronization of the SolidMechanicsModel.mass _gst_smm_for_strain, /// synchronization of the SolidMechanicsModel.current_position _gst_smm_boundary, /// synchronization of the boundary, forces, velocities and displacement /// Test tag _gst_test }; /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /// @enum SynchronizerOperation reduce operation that the synchronizer can perform enum SynchronizerOperation { _so_sum, _so_min, _so_max, _so_null }; /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name (type variable) { \ this->variable = variable; \ } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name () const { \ return variable; \ } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name () { \ return variable; \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ inline type get##name (const ::akantu::ElementType & el_type) const { \ AKANTU_DEBUG_IN(); \ AKANTU_DEBUG_ASSERT(variable[el_type] != NULL, \ "No " << #variable << " of type " \ << el_type << " in " << id); \ AKANTU_DEBUG_OUT(); \ return *variable[el_type]; \ } /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type) { switch(type) { case _segment_2 : stream << "segment_2" ; break; case _segment_3 : stream << "segment_3" ; break; case _triangle_3 : stream << "triangle_3" ; break; case _triangle_6 : stream << "triangle_6" ; break; case _tetrahedron_4 : stream << "tetrahedron_4" ; break; case _tetrahedron_10 : stream << "tetrahedron_10"; break; case _quadrangle_4 : stream << "quadrangle_4" ; break; + case _quadrangle_8 : stream << "quadrangle_8" ; break; case _hexahedron_8 : stream << "hexahedron_8" ; break; case _not_defined : stream << "undefined" ; break; case _max_element_type : stream << "ElementType(" << (int) type << ")"; break; case _point : stream << "point"; break; } return stream; } /* -------------------------------------------------------------------------- */ void initialize(int * argc, char *** argv); /* -------------------------------------------------------------------------- */ void finalize (); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ #if defined(__INTEL_COMPILER) /// remark #981: operands are evaluated in unspecified order #pragma warning ( disable : 981 ) /// remark #383: value copied to temporary, reference to temporary used #pragma warning ( disable : 383 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline void to_lower(std::string & str) { std::transform(str.begin(), str.end(), str.begin(), (int(*)(int))std::tolower); } /* -------------------------------------------------------------------------- */ inline void trim(std::string & to_trim) { size_t first = to_trim.find_first_not_of(" \t"); if (first != std::string::npos) { size_t last = to_trim.find_last_not_of(" \t"); to_trim = to_trim.substr(first, last - first + 1); } else to_trim = ""; } __END_AKANTU__ //#include "aka_types.hh" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/fem/element_class.cc b/fem/element_class.cc index 12e284c96..d9ccdc7de 100644 --- a/fem/element_class.cc +++ b/fem/element_class.cc @@ -1,203 +1,232 @@ /** * @file element_class.cc * @author Nicolas Richart * @date Tue Jul 20 10:12:44 2010 * * @brief Common part of element_classes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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<> UInt ElementClass<_point>::nb_quadrature_points = 1; template<> Real ElementClass<_point>::quad[] = {0}; 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<> 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[] = {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, 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[6], &vec_facet_connectivity[12], &vec_facet_connectivity[18]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_quadrangle_4>::nb_nodes_per_element = 4; template<> ElementType ElementClass<_quadrangle_4>::p1_element_type = _quadrangle_4; -template<> UInt ElementClass<_quadrangle_4>::nb_quadrature_points = 1; -template<> Real ElementClass<_quadrangle_4>::quad[] = {0, 0}; +template<> UInt ElementClass<_quadrangle_4>::nb_quadrature_points = 4; +//template<> Real ElementClass<_quadrangle_4>::quad[] = {0, 0}; +template<> Real ElementClass<_quadrangle_4>::quad[] = {-1./sqrt(3), -1./sqrt(3), + 1./sqrt(3), -1./sqrt(3), + 1./sqrt(3), 1./sqrt(3), + -1./sqrt(3), 1./sqrt(3)}; template<> UInt ElementClass<_quadrangle_4>::spatial_dimension = 2; template<> UInt ElementClass<_quadrangle_4>::nb_facets = 4; template<> ElementType ElementClass<_quadrangle_4>::facet_type = _segment_2; template<> UInt ElementClass<_quadrangle_4>::vec_facet_connectivity[]= {0, 1, 1, 2, 2, 3, 3, 0}; template<> UInt * ElementClass<_quadrangle_4>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[2], &vec_facet_connectivity[4], &vec_facet_connectivity[6]}; /* -------------------------------------------------------------------------- */ -// template<> UInt ElementClass<_quadrangle_6>::nb_nodes_per_element = 6; -// template<> ElementType ElementClass<_quadrangle_6>::p1_element_type = _quadrangle_6; -// template<> UInt ElementClass<_quadrangle_6>::nb_quadrature_points = 1; -// template<> Real ElementClass<_quadrangle_4>::quad[] = {0, 0}; -// template<> UInt ElementClass<_quadrangle_4>::spatial_dimension = 2; -// template<> UInt ElementClass<_quadrangle_4>::nb_facets = 4; -// template<> ElementType ElementClass<_quadrangle_4>::facet_type = _segment_2; -// template<> UInt ElementClass<_quadrangle_4>::vec_facet_connectivity[]= {0, 1, -// 1, 2, -// 2, 3, -// 3, 0}; -// template<> UInt * ElementClass<_quadrangle_4>::facet_connectivity[] = {&vec_facet_connectivity[0], -// &vec_facet_connectivity[2], -// &vec_facet_connectivity[4], -// &vec_facet_connectivity[6]}; +template<> UInt ElementClass<_quadrangle_8>::nb_nodes_per_element = 6; +template<> ElementType ElementClass<_quadrangle_8>::p1_element_type = _quadrangle_8; +template<> UInt ElementClass<_quadrangle_8>::nb_quadrature_points = 9; +template<> Real ElementClass<_quadrangle_8>::quad[] = { 0., 0., + sqrt(3./5.), sqrt(3./5.), + -sqrt(3./5.), sqrt(3./5.), + -sqrt(3./5.), -sqrt(3./5.), + sqrt(3./5.), -sqrt(3./5.), + 0., sqrt(3./5.), + -sqrt(3./5.), 0., + 0., -sqrt(3./5.), + sqrt(3./5.), 0.}; +template<> UInt ElementClass<_quadrangle_8>::spatial_dimension = 2; +template<> UInt ElementClass<_quadrangle_8>::nb_facets = 4; +template<> ElementType ElementClass<_quadrangle_8>::facet_type = _segment_2; +template<> UInt ElementClass<_quadrangle_8>::vec_facet_connectivity[]= {0, 1, 4, + 1, 2, 5, + 2, 3, 6, + 3, 0, 7}; +template<> UInt * ElementClass<_quadrangle_8>::facet_connectivity[] = {vec_facet_connectivity + 0, + vec_facet_connectivity + 3, + vec_facet_connectivity + 6, + vec_facet_connectivity + 9}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_hexahedron_8>::nb_nodes_per_element = 8; template<> ElementType ElementClass<_hexahedron_8>::p1_element_type = _hexahedron_8; template<> UInt ElementClass<_hexahedron_8>::nb_quadrature_points = 8; template<> Real ElementClass<_hexahedron_8>::quad[] = {-1./sqrt(3.), -1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.), -1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.), 1./sqrt(3.), -1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.), -1./sqrt(3.), -1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.), 1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.), 1./sqrt(3.), 1./sqrt(3.), 1./sqrt(3.), -1./sqrt(3.), 1./sqrt(3.), 1./sqrt(3.)}; template<> UInt ElementClass<_hexahedron_8>::spatial_dimension = 3; template<> UInt ElementClass<_hexahedron_8>::nb_facets = 6; template<> ElementType ElementClass<_hexahedron_8>::facet_type = _quadrangle_4; template<> UInt ElementClass<_hexahedron_8>::vec_facet_connectivity[]= {0, 1, 2, 3, 0, 1, 5, 4, 1, 2, 6, 5, 2, 3, 7, 6, 3, 0, 4, 7, 4, 5, 6, 7}; template<> UInt * ElementClass<_hexahedron_8>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[4], &vec_facet_connectivity[8], &vec_facet_connectivity[12], &vec_facet_connectivity[16], &vec_facet_connectivity[20]}; /* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- */ +/* Gauss integration */ +/* -------------------------------------------------------------------------- */ +template<> Real ElementClass<_segment_2 >::gauss_integration_weights[] = {2.}; +template<> Real ElementClass<_segment_3 >::gauss_integration_weights[] = {1., 1.}; +template<> Real ElementClass<_triangle_3 >::gauss_integration_weights[] = {1./2.}; +template<> Real ElementClass<_triangle_6 >::gauss_integration_weights[] = {1./6., 1./6., 1./6.}; +template<> Real ElementClass<_tetrahedron_4 >::gauss_integration_weights[] = {1./6.}; +template<> Real ElementClass<_tetrahedron_10>::gauss_integration_weights[] = {1./24., 1./24., 1./24., 1./24.}; +template<> Real ElementClass<_quadrangle_4 >::gauss_integration_weights[] = {1., 1., 1., 1.}; +template<> Real ElementClass<_quadrangle_8 >::gauss_integration_weights[] = {64./80., + 25./81., 25./81., 25./81., 25./81., + 40./81., 40./81., 40./81., 40./81.}; +template<> Real ElementClass<_hexahedron_8 >::gauss_integration_weights[] = {1., 1., 1., 1., + 1., 1., 1., 1.}; +template<> Real ElementClass<_point >::gauss_integration_weights[] = {1.}; + __END_AKANTU__ diff --git a/fem/element_class.hh b/fem/element_class.hh index 8891a6ea9..8a384ab5b 100644 --- a/fem/element_class.hh +++ b/fem/element_class.hh @@ -1,198 +1,203 @@ /** * @file element_class.hh * @author Nicolas Richart * @date Wed Jun 16 18:48:25 2010 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_HH__ #define __AKANTU_ELEMENT_CLASS_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Class describing the different type of element for mesh or finite element * purpose * * @tparam type the element type for the specialization of the element class */ template class ElementClass { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /** * compute the shape functions, the shape functions derivatives and the * jacobians * @param[in] coord coordinates of the nodes * @param[out] shape shape functions [nb_quad*node_per_elem] * @param[out] shape_deriv shape functions derivatives [nb_quad*node_per_elem*spatial_dim] * @param[out] jacobian jacobians * integration weights [nb_quad] */ inline static void preComputeStandards(const Real * coord, const UInt dimension, Real * shape, Real * shape_deriv, Real * jacobian); /// compute the shape values for a point given in natural coordinates inline static void computeShapes(const Real * natural_coords, Real * shapes); /// compute the shape values for a set of points given in natural coordinates inline static void computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes); /** compute dxds the variation of real coordinates along with variation of natural coordinates * on a given point in natural coordinates */ inline static void computeDXDS(const Real * dnds, const Real * node_coords, const UInt dimension, Real * dxds); /** compute dxds the variation of real coordinates along with variation of natural coordinates * on a given set of points in natural coordinates */ inline static void computeDXDS(const Real * dnds, const UInt nb_points, const Real * node_coords, const UInt dimension, Real * dxds); /** compute dnds the variation of real shape functions along with variation of natural coordinates * on a given point in natural coordinates */ inline static void computeDNDS(const Real * natural_coords, Real * dnds); /** compute dnds the variation of shape functions along with variation of natural coordinates * on a given set of points in natural coordinates */ inline static void computeDNDS(const Real * natural_coords, const UInt nb_points, Real * dnds); /// compute jacobian (or integration variable change factor) for a set of points inline static void computeJacobian(const Real * dxds, const UInt nb_points, const UInt dimension, Real * jac); /// compute jacobian (or integration variable change factor) for a given point inline static void computeJacobian(const Real * dxds, const UInt dimension, Real & jac); /// compute shape derivatives (input is dxds) for a set of points inline static void computeShapeDerivatives(const Real * dxds, const Real * dnds, const UInt nb_points, const UInt dimension, Real * shape_deriv); /// compute shape derivatives (input is dxds) for a given point inline static void computeShapeDerivatives(const Real * dxds, const Real * dnds, Real * shape_deriv); /// compute normals on quad points inline static void computeNormalsOnQuadPoint(const Real * dxds, const UInt dimension, Real * normals); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const {}; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: static AKANTU_GET_MACRO_NOT_CONST(NbNodesPerElement, nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, p1_element_type, ElementType); static AKANTU_GET_MACRO_NOT_CONST(NbQuadraturePoints, nb_quadrature_points, UInt); static AKANTU_GET_MACRO_NOT_CONST(SpatialDimension, spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(FacetElementType, facet_type, const ElementType &); static AKANTU_GET_MACRO_NOT_CONST(NbFacetsPerElement, nb_facets, UInt); static AKANTU_GET_MACRO_NOT_CONST(FacetLocalConnectivityPerElement, facet_connectivity, UInt**); static inline UInt getNbQuadraturePoint(); static inline Real * getQuadraturePoints(); static inline UInt getShapeSize(); static inline UInt getShapeDerivativesSize(); /// compute the in-radius static inline Real getInradius(const Real * coord); + static inline Real * getGaussIntegrationWeights(); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Number of nodes per element static UInt nb_nodes_per_element; /// Number of quadrature points per element static UInt nb_quadrature_points; /// Dimension of the element static UInt spatial_dimension; /// Type of the facet elements static ElementType facet_type; /// number of facets for element static UInt nb_facets; /// local connectivity of facets static UInt * facet_connectivity[]; /// vectorial connectivity of facets static UInt vec_facet_connectivity[]; /// type of element P1 associated static ElementType p1_element_type; /// quadrature points in natural coordinates static Real quad[]; + + /// weights for the Gauss integration + static Real gauss_integration_weights[]; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "element_class_inline_impl.cc" /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_ELEMENT_CLASS_HH__ */ diff --git a/fem/element_class_inline_impl.cc b/fem/element_class_inline_impl.cc index df1bb55b0..294425cca 100644 --- a/fem/element_class_inline_impl.cc +++ b/fem/element_class_inline_impl.cc @@ -1,244 +1,249 @@ /** * @file element_class_inline_impl.cc * @author Nicolas Richart * @author Guillaume Anciaux * @date Thu Jul 15 10:28:28 2010 * * @brief Implementation of the inline functions of the class element_class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getNbQuadraturePoint() { return nb_quadrature_points; } /* -------------------------------------------------------------------------- */ template inline Real * ElementClass::getQuadraturePoints() { return quad; } /* -------------------------------------------------------------------------- */ 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 + // 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, +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, + 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 +template inline void ElementClass::computeDXDS(const Real * dnds, - const Real * node_coords, - const UInt dimension, + 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 +template inline void ElementClass::computeJacobian(const Real * dxds, const UInt nb_points, - const UInt dimension, + 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, + const UInt dimension, Real * shape_deriv) { - AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"gradient in space " - << dimension + 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$ 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); - + "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); + 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){ +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 +template inline void ElementClass::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, - __attribute__ ((unused)) Real * dnds){ + __attribute__ ((unused)) Real * dnds) { AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ -template +template inline void ElementClass::computeJacobian(__attribute__ ((unused)) const Real * dxds, - __attribute__ ((unused)) const UInt dimension, - __attribute__ ((unused)) Real & jac){ + __attribute__ ((unused)) const UInt dimension, + __attribute__ ((unused)) Real & jac) { AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } +/* -------------------------------------------------------------------------- */ +template +inline Real * ElementClass::getGaussIntegrationWeights() { + return gauss_integration_weights; +} #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_quadrangle_4.cc" #include "element_classes/element_class_hexahedron_8.cc" diff --git a/fem/element_classes/element_class_hexahedron_8.cc b/fem/element_classes/element_class_hexahedron_8.cc index 685999506..b859e5e54 100644 --- a/fem/element_classes/element_class_hexahedron_8.cc +++ b/fem/element_classes/element_class_hexahedron_8.cc @@ -1,211 +1,210 @@ /** * @file element_class_hexahedron_8.cc * @author Peter Spijker * @date Mon Mar 14 17:27:44 2010 * * @brief Specialization of the element_class class for the type _hexahedron_8 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @verbatim \zeta ^ (-1,1,1) | (1,1,1) 8---|------7 /| | /| / | | / | (-1,-1,1) 5----------6 | (1,-1,1) | | | | | | | | | | | | +---|-------> \xi | | / | | (-1,1,-1) | 4-/-----|--3 (1,1,-1) | / / | / |/ / |/ 1-/--------2 (-1,-1,-1) / (1,-1,-1) / \eta @endverbatim * * @subsection shapes Shape functions * @f[ * \begin{array}{llll} * N1 = (1 - \xi) (1 - \eta) (1 - \zeta) / 8 * & \frac{\partial N1}{\partial \xi} = - (1 - \eta) (1 - \zeta) / 8 * & \frac{\partial N1}{\partial \eta} = - (1 - \xi) (1 - \zeta) / 8 * & \frac{\partial N1}{\partial \zeta} = - (1 - \xi) (1 - \eta) / 8 \\ * N2 = (1 + \xi) (1 - \eta) (1 - \zeta) / 8 * & \frac{\partial N2}{\partial \xi} = (1 - \eta) (1 - \zeta) / 8 * & \frac{\partial N2}{\partial \eta} = - (1 + \xi) (1 - \zeta) / 8 * & \frac{\partial N2}{\partial \zeta} = - (1 + \xi) (1 - \eta) / 8 \\ * N3 = (1 + \xi) (1 + \eta) (1 - \zeta) / 8 * & \frac{\partial N3}{\partial \xi} = (1 + \eta) (1 - \zeta) / 8 * & \frac{\partial N3}{\partial \eta} = (1 + \xi) (1 - \zeta) / 8 * & \frac{\partial N3}{\partial \zeta} = - (1 + \xi) (1 + \eta) / 8 \\ * N4 = (1 - \xi) (1 + \eta) (1 - \zeta) / 8 * & \frac{\partial N4}{\partial \xi} = - (1 + \eta) (1 - \zeta) / 8 * & \frac{\partial N4}{\partial \eta} = (1 - \xi) (1 - \zeta) / 8 * & \frac{\partial N4}{\partial \zeta} = - (1 - \xi) (1 + \eta) / 8 \\ * N5 = (1 - \xi) (1 - \eta) (1 + \zeta) / 8 * & \frac{\partial N5}{\partial \xi} = - (1 - \eta) (1 + \zeta) / 8 * & \frac{\partial N5}{\partial \eta} = - (1 - \xi) (1 + \zeta) / 8 * & \frac{\partial N5}{\partial \zeta} = (1 - \xi) (1 - \eta) / 8 \\ * N6 = (1 + \xi) (1 - \eta) (1 + \zeta) / 8 * & \frac{\partial N6}{\partial \xi} = (1 - \eta) (1 + \zeta) / 8 * & \frac{\partial N6}{\partial \eta} = - (1 + \xi) (1 + \zeta) / 8 * & \frac{\partial N6}{\partial \zeta} = (1 + \xi) (1 - \eta) / 8 \\ * N7 = (1 + \xi) (1 + \eta) (1 + \zeta) / 8 * & \frac{\partial N7}{\partial \xi} = (1 + \eta) (1 + \zeta) / 8 * & \frac{\partial N7}{\partial \eta} = (1 + \xi) (1 + \zeta) / 8 * & \frac{\partial N7}{\partial \zeta} = (1 + \xi) (1 + \eta) / 8 \\ * N8 = (1 - \xi) (1 + \eta) (1 + \zeta) / 8 * & \frac{\partial N8}{\partial \xi} = - (1 + \eta) (1 + \zeta) / 8 * & \frac{\partial N8}{\partial \eta} = (1 - \xi) (1 + \zeta) / 8 * & \frac{\partial N8}{\partial \zeta} = (1 - \xi) (1 + \eta) / 8 \\ * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& -1/\sqrt{3} \qquad \eta_{q0} = -1/\sqrt{3} \qquad \zeta_{q0} = -1/\sqrt{3} \\ * \xi_{q1} &=& 1/\sqrt{3} \qquad \eta_{q1} = -1/\sqrt{3} \qquad \zeta_{q1} = -1/\sqrt{3} \\ * \xi_{q2} &=& 1/\sqrt{3} \qquad \eta_{q2} = 1/\sqrt{3} \qquad \zeta_{q2} = -1/\sqrt{3} \\ * \xi_{q3} &=& -1/\sqrt{3} \qquad \eta_{q3} = 1/\sqrt{3} \qquad \zeta_{q3} = -1/\sqrt{3} \\ * \xi_{q4} &=& -1/\sqrt{3} \qquad \eta_{q4} = -1/\sqrt{3} \qquad \zeta_{q4} = 1/\sqrt{3} \\ * \xi_{q5} &=& 1/\sqrt{3} \qquad \eta_{q5} = -1/\sqrt{3} \qquad \zeta_{q5} = 1/\sqrt{3} \\ * \xi_{q6} &=& 1/\sqrt{3} \qquad \eta_{q6} = 1/\sqrt{3} \qquad \zeta_{q6} = 1/\sqrt{3} \\ * \xi_{q7} &=& -1/\sqrt{3} \qquad \eta_{q7} = 1/\sqrt{3} \qquad \zeta_{q7} = 1/\sqrt{3} \\ * @f} */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_hexahedron_8>::nb_nodes_per_element; template<> UInt ElementClass<_hexahedron_8>::nb_quadrature_points; template<> UInt ElementClass<_hexahedron_8>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_hexahedron_8>::computeShapes(const Real * natural_coords, Real * shapes) { /// Natural coordinates const Real * c = natural_coords; shapes[0] = .125 * (1 - c[0]) * (1 - c[1]) * (1 - c[3]); /// N1(q_0) shapes[1] = .125 * (1 + c[0]) * (1 - c[1]) * (1 - c[3]); /// N2(q_0) shapes[2] = .125 * (1 + c[0]) * (1 + c[1]) * (1 - c[3]); /// N3(q_0) shapes[3] = .125 * (1 - c[0]) * (1 + c[1]) * (1 - c[3]); /// N4(q_0) shapes[4] = .125 * (1 - c[0]) * (1 - c[1]) * (1 + c[3]); /// N5(q_0) shapes[5] = .125 * (1 + c[0]) * (1 - c[1]) * (1 + c[3]); /// N6(q_0) shapes[6] = .125 * (1 + c[0]) * (1 + c[1]) * (1 + c[3]); /// N7(q_0) shapes[7] = .125 * (1 - c[0]) * (1 + c[1]) * (1 + c[3]); /// N8(q_0) } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_hexahedron_8>::computeDNDS(const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccccccc} * \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 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 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} * \end{array} * \right) * @f] */ const Real * c = natural_coords; dnds[0] = - .125 * (1 - c[1]) * (1 - c[2]);; dnds[1] = .125 * (1 - c[1]) * (1 - c[2]);; dnds[2] = .125 * (1 + c[1]) * (1 - c[2]);; dnds[3] = - .125 * (1 + c[1]) * (1 - c[2]);; dnds[4] = - .125 * (1 - c[1]) * (1 + c[2]);; dnds[5] = .125 * (1 - c[1]) * (1 + c[2]);; dnds[6] = .125 * (1 + c[1]) * (1 + c[2]);; dnds[7] = - .125 * (1 + c[1]) * (1 + c[2]);; dnds[8] = - .125 * (1 - c[0]) * (1 - c[2]);; dnds[9] = - .125 * (1 + c[0]) * (1 - c[2]);; dnds[10] = .125 * (1 + c[0]) * (1 - c[2]);; dnds[11] = .125 * (1 - c[0]) * (1 - c[2]);; dnds[12] = - .125 * (1 - c[0]) * (1 + c[2]);; dnds[13] = - .125 * (1 + c[0]) * (1 + c[2]);; dnds[14] = .125 * (1 + c[0]) * (1 + c[2]);; dnds[15] = .125 * (1 - c[0]) * (1 + c[2]);; dnds[16] = - .125 * (1 - c[0]) * (1 - c[1]);; dnds[17] = - .125 * (1 + c[0]) * (1 - c[1]);; dnds[18] = - .125 * (1 + c[0]) * (1 + c[1]);; dnds[19] = - .125 * (1 - c[0]) * (1 + c[1]);; dnds[20] = .125 * (1 - c[0]) * (1 - c[1]);; dnds[21] = .125 * (1 + c[0]) * (1 - c[1]);; dnds[22] = .125 * (1 + c[0]) * (1 + c[1]);; dnds[23] = .125 * (1 - c[0]) * (1 + c[1]);; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_hexahedron_8>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ - Real weight = 1.; if (dimension == spatial_dimension){ Real det_dxds = Math::det3(dxds); - jac = det_dxds * weight; + jac = det_dxds; } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_hexahedron_8>::getInradius(const Real * coord) { Real a = Math::distance_3d(coord + 0, coord + 3); Real b = Math::distance_3d(coord + 3, coord + 6); Real c = Math::distance_3d(coord + 6, coord + 9); Real d = Math::distance_3d(coord + 9, coord + 0); Real e = Math::distance_3d(coord + 0, coord + 12); Real f = Math::distance_3d(coord + 3, coord + 15); Real g = Math::distance_3d(coord + 6, coord + 18); Real h = Math::distance_3d(coord + 9, coord + 21); Real i = Math::distance_3d(coord + 12, coord + 15); Real j = Math::distance_3d(coord + 15, coord + 18); Real k = Math::distance_3d(coord + 18, coord + 21); Real l = Math::distance_3d(coord + 21, coord + 12); Real x = std::min(a, std::min(b, std::min(c, d))); Real y = std::min(e, std::min(f, std::min(g, h))); Real z = std::min(i, std::min(j, std::min(k, l))); Real p = std::min(x, std::min(y, z)); return p; } diff --git a/fem/element_classes/element_class_quadrangle_4.cc b/fem/element_classes/element_class_quadrangle_4.cc index 7b3aa683f..21a636614 100644 --- a/fem/element_classes/element_class_quadrangle_4.cc +++ b/fem/element_classes/element_class_quadrangle_4.cc @@ -1,145 +1,144 @@ /** * @file element_class_quadrangle_4.cc * @author Nicolas Richart * @date Wed Oct 27 17:27:44 2010 * * @brief Specialization of the element_class class for the type _quadrangle_4 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @verbatim \eta ^ (-1,1) | (1,1) x---------x | | | | | | --|---------|-----> \xi | | | | | | x---------x (-1,-1) | (1,-1) @endverbatim * * @subsection shapes Shape functions * @f[ * \begin{array}{lll} * N1 = (1 - \xi) (1 - \eta) / 4 * & \frac{\partial N1}{\partial \xi} = - (1 - \eta) / 4 * & \frac{\partial N1}{\partial \eta} = - (1 - \xi) / 4 \\ * N2 = (1 + \xi) (1 - \eta) / 4 \\ * & \frac{\partial N2}{\partial \xi} = (1 - \eta) / 4 * & \frac{\partial N2}{\partial \eta} = - (1 + \xi) / 4 \\ * N3 = (1 + \xi) (1 + \eta) / 4 \\ * & \frac{\partial N3}{\partial \xi} = (1 + \eta) / 4 * & \frac{\partial N3}{\partial \eta} = (1 + \xi) / 4 \\ * N4 = (1 - \xi) (1 + \eta) / 4 * & \frac{\partial N4}{\partial \xi} = - (1 + \eta) / 4 * & \frac{\partial N4}{\partial \eta} = (1 - \xi) / 4 \\ * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& 0 \qquad \eta_{q0} = 0 * @f} */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_quadrangle_4>::nb_nodes_per_element; template<> UInt ElementClass<_quadrangle_4>::nb_quadrature_points; template<> UInt ElementClass<_quadrangle_4>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeShapes(const Real * natural_coords, Real * shapes) { /// Natural coordinates const Real * c = natural_coords; shapes[0] = .25 * (1 - c[0]) * (1 - c[1]); /// N1(q_0) shapes[1] = .25 * (1 + c[0]) * (1 - c[1]); /// N2(q_0) shapes[2] = .25 * (1 + c[0]) * (1 + c[1]); /// N3(q_0) shapes[3] = .25 * (1 - c[0]) * (1 + c[1]); /// N4(q_0) } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeDNDS(const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccc} * \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} * \end{array} * \right) * @f] */ const Real * c = natural_coords; dnds[0] = - .25 * (1 - c[1]); dnds[1] = .25 * (1 - c[1]); dnds[2] = .25 * (1 + c[1]); dnds[3] = - .25 * (1 + c[1]); dnds[4] = - .25 * (1 - c[0]); dnds[5] = - .25 * (1 + c[0]); dnds[6] = .25 * (1 + c[0]); dnds[7] = .25 * (1 - c[0]); } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ - Real weight = 4.; if (dimension == spatial_dimension){ Real det_dxds = Math::det2(dxds); - jac = det_dxds * weight; + jac = det_dxds; } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_quadrangle_4>::getInradius(const Real * coord) { Real a = Math::distance_2d(coord + 0, coord + 2); Real b = Math::distance_2d(coord + 2, coord + 4); Real c = Math::distance_2d(coord + 4, coord + 6); Real d = Math::distance_2d(coord + 6, coord + 0); // Real septimetre = (a + b + c + d) / 2.; // Real p = Math::distance_2d(coord + 0, coord + 4); // Real q = Math::distance_2d(coord + 2, coord + 6); // Real area = sqrt(4*(p*p * q*q) - (a*a + b*b + c*c + d*d)*(a*a + c*c - b*b - d*d)) / 4.; // Real h = sqrt(area); // to get a length // Real h = area / septimetre; // formula of inradius for circumscritable quadrelateral Real h = std::min(a, std::min(b, std::min(c, d))); return h; } diff --git a/fem/element_classes/element_class_quadrangle_8.cc b/fem/element_classes/element_class_quadrangle_8.cc new file mode 100644 index 000000000..f86d05f44 --- /dev/null +++ b/fem/element_classes/element_class_quadrangle_8.cc @@ -0,0 +1,151 @@ +/** + * @file element_class_quadrangle_8.cc + * @author Nicolas Richart + * @date Tue May 17 18:35:04 2011 + * + * @brief Specialization of the ElementClass for the _quadrangle_8 + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + * @verbatim + \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{equerry*}{ + * \xi_{q0} &=& 0 \qquad \eta_{q0} = 0 + * @f} + */ + +/* -------------------------------------------------------------------------- */ +template<> UInt ElementClass<_quadrangle_8>::nb_nodes_per_element; +template<> UInt ElementClass<_quadrangle_8>::nb_quadrature_points; +template<> UInt ElementClass<_quadrangle_8>::spatial_dimension; + +/* -------------------------------------------------------------------------- */ +template <> inline void ElementClass<_quadrangle_8>::computeShapes(const Real * natural_coords, + Real * shapes) { + /// Natural coordinates + const Real * c = natural_coords; + + shapes[0] = - .25 * (1 - c[0]) * (1 - c[1]) * (1 + c[0] + c[1]); + shapes[1] = - .25 * (1 + c[0]) * (1 - c[1]) * (1 - c[0] + c[1]); + shapes[2] = - .25 * (1 + c[0]) * (1 + c[1]) * (1 - c[0] - c[1]); + shapes[3] = - .25 * (1 - c[0]) * (1 + c[1]) * (1 + c[0] - c[1]); + shapes[4] = .5 * (1 - c[0] * c[0]) * (1 - c[1] ); + shapes[5] = .5 * (1 + c[0] ) * (1 - c[1] * c[1]); + shapes[6] = .5 * (1 - c[0] * c[0]) * (1 + c[1] ); + shapes[7] = .5 * (1 - c[0] ) * (1 - c[1] * c[1]); +} +/* -------------------------------------------------------------------------- */ +template <> inline void ElementClass<_quadrangle_8>::computeDNDS(const Real * natural_coords, + Real * dnds) { + const Real * c = natural_coords; + + dnds[0] = .25 * (1 - c[1]) * (2 * c[0] + c[1]); + dnds[1] = .25 * (1 - c[1]) * (2 * c[0] - c[1]); + dnds[2] = .25 * (1 + c[1]) * (2 * c[0] + c[1]); + dnds[3] = .25 * (1 + c[1]) * (2 * c[0] - c[1]); + dnds[4] = - c[0] * (1 - c[1]); + dnds[5] = .5 * (1 - c[1] * c[1]); + dnds[6] = - c[0] * (1 + c[1]); + dnds[7] = - .5 * (1 - c[1] * c[1]); + + dnds[8] = .25 * (1 - c[0]) * (2 * c[1] + c[0]); + dnds[9] = .25 * (1 + c[0]) * (2 * c[1] - c[0]); + dnds[10] = .25 * (1 + c[0]) * (2 * c[1] + c[0]); + dnds[11] = .25 * (1 - c[0]) * (2 * c[1] - c[0]); + dnds[12] = - .5 * (1 - c[0] * c[0]); + dnds[13] = - c[1] * (1 + c[0]); + dnds[14] = .5 * (1 - c[0] * c[0]); + dnds[15] = - c[1] * (1 - c[0]); +} + + +/* -------------------------------------------------------------------------- */ +template <> inline void ElementClass<_quadrangle_8>::computeJacobian(const Real * dxds, + const UInt dimension, + Real & jac){ + if (dimension == spatial_dimension){ + Real det_dxds = Math::det2(dxds); + jac = det_dxds; + } else { + AKANTU_DEBUG_TO_IMPLEMENT(); + } +} + + + +/* -------------------------------------------------------------------------- */ +template<> inline Real ElementClass<_quadrangle_4>::getInradius(const Real * coord) { + Real a = Math::distance_2d(coord + 0, coord + 2); + Real b = Math::distance_2d(coord + 2, coord + 4); + Real c = Math::distance_2d(coord + 4, coord + 6); + Real d = Math::distance_2d(coord + 6, coord + 0); + + Real h = std::min(a, std::min(b, std::min(c, d))); + + return h; +} diff --git a/fem/element_classes/element_class_segment_2.cc b/fem/element_classes/element_class_segment_2.cc index 392c8f628..18e8067c2 100644 --- a/fem/element_classes/element_class_segment_2.cc +++ b/fem/element_classes/element_class_segment_2.cc @@ -1,92 +1,90 @@ /** * @file element_class_segment_2.cc * @author Nicolas Richart * @date Thu Jul 15 10:28:28 2010 * * @brief Specialization of the element_class class for the type _segment_2 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @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} */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_segment_2>::nb_nodes_per_element; template<> UInt ElementClass<_segment_2>::nb_quadrature_points; template<> UInt ElementClass<_segment_2>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_segment_2>::computeShapes(const Real * natural_coords, Real * shapes){ /// natural coordinate Real c = natural_coords[0]; /// shape functions shapes[0] = 0.5*(1-c); shapes[1] =0.5*(1+c); } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_segment_2>::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, Real * dnds){ /// dN1/de dnds[0] = - .5; /// dN2/de dnds[1] = .5; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_segment_2>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ if (dimension == spatial_dimension){ jac = dxds[0]; } else { jac = Math::norm2(dxds); } - // multiplication by integration weight - jac *= 2; } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_segment_2>::getInradius(const Real * coord) { return sqrt((coord[0] - coord[1])*(coord[0] - coord[1])); } diff --git a/fem/element_classes/element_class_segment_3.cc b/fem/element_classes/element_class_segment_3.cc index d4b9347a0..62fa01c76 100644 --- a/fem/element_classes/element_class_segment_3.cc +++ b/fem/element_classes/element_class_segment_3.cc @@ -1,107 +1,103 @@ /** * @file element_class_segment_3.cc * @author Nicolas Richart * @date Sun Oct 3 10:28:28 2010 * * @brief Specialization of the element_class class for the type _segment_3 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @verbatim -1 0 1 -----x---------x---------x-----> x 1 3 2 @endverbatim * * @subsection coords Nodes coordinates * * @f[ * \begin{array}{lll} * x_{1} = -1 & x_{2} = 1 & x_{3} = 0 * \end{array} * @f] * * @subsection shapes Shape functions * @f[ * \begin{array}{ll} * w_1(x) = \frac{x}{2}(x - 1) & w'_1(x) = x - \frac{1}{2}\\ * w_2(x) = \frac{x}{2}(x + 1) & w'_2(x) = x + \frac{1}{2}\\ * w_3(x) = 1-x^2 & w'_3(x) = -2x * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f[ * \begin{array}{ll} * x_{q1} = -1/\sqrt{3} & x_{q2} = 1/\sqrt{3} * \end{array} * @f] */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_segment_3>::nb_nodes_per_element; template<> UInt ElementClass<_segment_3>::nb_quadrature_points; template<> UInt ElementClass<_segment_3>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_segment_3>::computeShapes(const Real * natural_coords, Real * shapes){ Real c = natural_coords[0]; shapes[0] = (c - 1) * c / 2; shapes[1] = (c + 1) * c / 2; shapes[2] = 1 - c * c; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_segment_3>::computeDNDS(const Real * natural_coords, Real * dnds){ Real c = natural_coords[0]; dnds[0] = c - .5; dnds[1] = c + .5; dnds[2] = -2 * c; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_segment_3>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ - Real weight = 1; if (dimension == spatial_dimension){ jac = dxds[0]; - } - else { + } else { jac = Math::norm2(dxds); } - - jac *= weight; } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_segment_3>::getInradius(const Real * coord) { Real dist1 = sqrt((coord[0] - coord[1])*(coord[0] - coord[1])); Real dist2 = sqrt((coord[1] - coord[2])*(coord[1] - coord[2])); return dist1 < dist2 ? dist1 : dist2; } diff --git a/fem/element_classes/element_class_tetrahedron_10.cc b/fem/element_classes/element_class_tetrahedron_10.cc index 9bae4634b..b12ed5a83 100644 --- a/fem/element_classes/element_class_tetrahedron_10.cc +++ b/fem/element_classes/element_class_tetrahedron_10.cc @@ -1,269 +1,268 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @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) - ` ' | 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} = 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} = 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} = 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} = 4 \zeta - 1 \\ * N5 = 4 \xi (1 - \xi - \eta - \zeta) * & \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 \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} = 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, +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 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, + const UInt dimension, Real & jac) { if (dimension == spatial_dimension){ - Real weight = 1./24.; Real det_dxds = Math::det3(dxds); - jac = det_dxds * weight; - } + jac = det_dxds; + } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } - + /* -------------------------------------------------------------------------- */ 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 430770c35..8af494698 100644 --- a/fem/element_classes/element_class_tetrahedron_4.cc +++ b/fem/element_classes/element_class_tetrahedron_4.cc @@ -1,132 +1,130 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @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; + jac = det_dxds; } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_tetrahedron_4>::getInradius(const Real * coord) { return Math::tetrahedron_inradius(coord, coord+3, coord+6, coord+9); } diff --git a/fem/element_classes/element_class_triangle_3.cc b/fem/element_classes/element_class_triangle_3.cc index cc5140058..75c6bad13 100644 --- a/fem/element_classes/element_class_triangle_3.cc +++ b/fem/element_classes/element_class_triangle_3.cc @@ -1,122 +1,120 @@ /** * @file element_class_triangle_3.cc * @author Nicolas Richart * @date Thu Jul 15 10:28:28 2010 * * @brief Specialization of the element_class class for the type _triangle_3 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @verbatim \eta ^ | x (0,0,1) |` | ` | q ` | ° ` x--------x-----> \xi (1,0,0) (0,1,0) @endverbatim * * @subsection shapes Shape functions * @f{eqnarray*}{ * N1 &=& 1 - \xi - \eta \\ * N2 &=& \xi \\ * N3 &=& \eta * @f} * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& 1/3 \qquad \eta_{q0} = 1/3 * @f} */ // /// shape functions // shape[0] = 1./3.; /// N1(q_0) // shape[1] = 1./3.; /// N2(q_0) // shape[2] = 1./3.; /// N3(q_0) /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_3>::nb_nodes_per_element; template<> UInt ElementClass<_triangle_3>::nb_quadrature_points; template<> UInt ElementClass<_triangle_3>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_3>::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; /// N1(q_0) shapes[1] = c1; /// N2(q_0) shapes[2] = c2; /// N3(q_0) } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_3>::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 N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} & \frac{\partial N3}{\partial \eta} * \end{array} * \right) * @f] */ dnds[0] = -1.; dnds[1] = 1.; dnds[2] = 0.; dnds[3] = -1.; dnds[4] = 0.; dnds[5] = 1.; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_3>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ - const Real weight = .5; if (dimension == spatial_dimension){ Real det_dxds = Math::det2(dxds); jac = det_dxds; } else { Real vprod[dimension]; Math::vectorProduct3(dxds,dxds+3,vprod); jac = Math::norm3(vprod); } - jac *= weight; } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_triangle_3>::getInradius(const Real * coord) { return Math::triangle_inradius(coord, coord+2, coord+4); } diff --git a/fem/element_classes/element_class_triangle_6.cc b/fem/element_classes/element_class_triangle_6.cc index ea17a1ebc..4ca72acb9 100644 --- a/fem/element_classes/element_class_triangle_6.cc +++ b/fem/element_classes/element_class_triangle_6.cc @@ -1,208 +1,207 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @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 N2}{\partial \xi} = - 1 + 4 \xi * & \frac{\partial N2}{\partial \eta} = 0 \\ * N3 = - \eta (1 - 2 \eta) * & \frac{\partial N3}{\partial \xi} = 0 * & \frac{\partial N3}{\partial \eta} = - 1 + 4 \eta \\ * N4 = 4 \xi (1 - \xi - \eta) * & \frac{\partial N4}{\partial \xi} = 4 (1 - 2 \xi - \eta) * & \frac{\partial N4}{\partial \eta} = - 4 \eta \\ * N5 = 4 \xi \eta * & \frac{\partial N5}{\partial \xi} = 4 \xi * & \frac{\partial N5}{\partial \eta} = 4 \eta \\ * N6 = 4 \eta (1 - \xi - \eta) * & \frac{\partial N6}{\partial \xi} = - 4 \xi * & \frac{\partial N6}{\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, +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) + * \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, + 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; + jac = Math::det2(dxds); 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] * spatial_dimension, coord + triangles[t][1] * spatial_dimension, coord + triangles[t][2] * spatial_dimension); inradius = ir < inradius ? ir : inradius; } return inradius; } diff --git a/fem/integrator_gauss.cc b/fem/integrator_gauss.cc index 86032c16c..9262230a0 100644 --- a/fem/integrator_gauss.cc +++ b/fem/integrator_gauss.cc @@ -1,253 +1,258 @@ /** * @file integrator_gauss.cc * @author Guillaume Anciaux * @author Nicolas Richart * @date Thu Feb 10 16:52:07 2011 * * @brief implementation of gauss integrator class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ IntegratorGauss::IntegratorGauss(Mesh & mesh, IntegratorID id) : Integrator(mesh,id){ for(UInt t = _not_defined; t < _max_element_type; ++t) { this->jacobians[t] = NULL; this->ghost_jacobians[t] = NULL; quadrature_points[t] = NULL; } } /* -------------------------------------------------------------------------- */ template void IntegratorGauss::precomputeJacobiansOnQuadraturePoints(const UInt dimension, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Real * coord = mesh->getNodes().values; UInt spatial_dimension = mesh->getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); + Real * weights = ElementClass::getGaussIntegrationWeights(); + UInt * elem_val; UInt nb_element; std::string ghost = ""; if(ghost_type == _not_ghost) { elem_val = mesh->getConnectivity(type).values; nb_element = mesh->getConnectivity(type).getSize(); } else { ghost = "ghost_"; elem_val = mesh->getGhostConnectivity(type).values; nb_element = mesh->getGhostConnectivity(type).getSize(); } std::stringstream sstr_jacobians; sstr_jacobians << id << ":" << ghost << "jacobians:" << type; Vector * jacobians_tmp = &(alloc(sstr_jacobians.str(), nb_element*nb_quadrature_points, 1)); Real * jacobians_val = jacobians_tmp->values; Real local_coord[spatial_dimension * nb_nodes_per_element]; for (UInt elem = 0; elem < nb_element; ++elem) { mesh->extractNodalCoordinatesFromElement(local_coord, elem_val+elem*nb_nodes_per_element, nb_nodes_per_element); Real * quad = ElementClass::getQuadraturePoints(); // compute dnds Real dnds[nb_nodes_per_element * spatial_dimension * nb_quadrature_points]; ElementClass::computeDNDS(quad, nb_quadrature_points, dnds); // compute dxds Real dxds[dimension * spatial_dimension * nb_quadrature_points]; ElementClass::computeDXDS(dnds, nb_quadrature_points, local_coord, dimension, dxds); // jacobian ElementClass::computeJacobian(dxds, nb_quadrature_points, dimension, jacobians_val); - jacobians_val += nb_quadrature_points; + for (UInt q = 0; q < nb_quadrature_points; ++q) { + *jacobians_val++ *= weights[q]; + } + // jacobians_val += nb_quadrature_points; } if(ghost_type == _not_ghost) jacobians[type] = jacobians_tmp; else ghost_jacobians[type] = jacobians_tmp; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void IntegratorGauss::integrate(const Vector & in_f, Vector &intf, UInt nb_degre_of_freedom, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * jac_loc; UInt nb_element; if(ghost_type == _not_ghost) { jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points, "The vector in_f(" << in_f.getID() << " size " << in_f.getSize() << ") has not the good size (" << nb_element << ")."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == nb_degre_of_freedom , "The vector in_f(" << in_f.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degre_of_freedom, "The vector intf(" << intf.getID() << ") has not the good number of component."); intf.resize(nb_element); Real * in_f_val = in_f.values; Real * intf_val = intf.values; Real * jac_val = jac_loc->values; UInt offset_in_f = in_f.getNbComponent()*nb_quadrature_points; UInt offset_intf = intf.getNbComponent(); Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * nb_quadrature_points; } integrate(in_f_val, jac, intf_val, nb_degre_of_freedom, nb_quadrature_points); in_f_val += offset_in_f; intf_val += offset_intf; if(filter_elements == NULL) { jac += nb_quadrature_points; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real IntegratorGauss::integrate(const Vector & in_f, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * jac_loc; UInt nb_element; if(ghost_type == _not_ghost) { jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points, "The vector in_f(" << in_f.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == 1, "The vector in_f(" << in_f.getID() << ") has not the good number of component."); Real intf = 0.; Real * in_f_val = in_f.values; Real * jac_val = jac_loc->values; UInt offset_in_f = in_f.getNbComponent() * nb_quadrature_points; Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * nb_quadrature_points; } Real el_intf = 0; integrate(in_f_val, jac, &el_intf, 1, nb_quadrature_points); intf += el_intf; in_f_val += offset_in_f; if(filter_elements == NULL) { jac += nb_quadrature_points; } } AKANTU_DEBUG_OUT(); return intf; } /* -------------------------------------------------------------------------- */ /* template instanciation */ /* -------------------------------------------------------------------------- */ #define INSTANCIATE_TEMPLATE_CLASS(type) \ template void IntegratorGauss::precomputeJacobiansOnQuadraturePoints(const UInt dimension, \ GhostType ghost_type); \ template void IntegratorGauss::integrate(const Vector & in_f, \ Vector &intf, \ UInt nb_degre_of_freedom, \ GhostType ghost_type, \ const Vector * filter_elements) const; \ template Real IntegratorGauss::integrate(const Vector & in_f, \ GhostType ghost_type, \ const Vector * filter_elements) const; AKANTU_BOOST_ELEMENT_LIST(INSTANCIATE_TEMPLATE_CLASS) #undef INSTANCIATE_TEMPLATE_CLASS __END_AKANTU__ diff --git a/fem/integrator_gauss.hh b/fem/integrator_gauss.hh index 56fc51512..69d871f09 100644 --- a/fem/integrator_gauss.hh +++ b/fem/integrator_gauss.hh @@ -1,122 +1,122 @@ /** * @file integrator_gauss.hh * @author Guillaume Anciaux * @author Nicolas Richart * @date Thu Feb 10 16:42:34 2011 * * @brief Gauss integration facilities * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "integrator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class IntegratorGauss : public Integrator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: IntegratorGauss(Mesh & mesh, IntegratorID id="IntegratorGauss"); virtual ~IntegratorGauss(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// precompute jacobians on elements of type "type" template void precomputeJacobiansOnQuadraturePoints(const UInt dimension, GhostType ghost_type); /// integrate f on the element "elem" of type "type" template inline void integrateOnElement(const Vector & f, Real * intf, UInt nb_degre_of_freedom, const UInt elem, GhostType ghost_type) const; /// integrate f for all elements of type "type" template void integrate(const Vector & in_f, Vector &intf, UInt nb_degre_of_freedom, GhostType ghost_type, const Vector * filter_elements) const; /// integrate scalar field in_f template Real integrate(const Vector & in_f, GhostType ghost_type, const Vector * filter_elements) const; /// return a vector with quadrature points natural coordinates template Vector & getQuadraturePoints() const; /// compute the vector of quadrature points natural coordinates template void computeQuadraturePoints(); // /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ - +private: // /// get the number of quadrature points // static inline UInt getNbQuadraturePoints(const ElementType & type); public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: inline void integrate(Real *f, Real *jac, Real * inte, UInt nb_degre_of_freedom, UInt nb_quadrature_points) const; ByElementTypeReal quadrature_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const IntegratorGauss & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ diff --git a/mesh_utils/mesh_io/mesh_io_msh.cc b/mesh_utils/mesh_io/mesh_io_msh.cc index b12c69000..5e7a14348 100644 --- a/mesh_utils/mesh_io/mesh_io_msh.cc +++ b/mesh_utils/mesh_io/mesh_io_msh.cc @@ -1,433 +1,435 @@ /** * @file mesh_io_msh.cc * @author Nicolas Richart * @date Fri Jun 18 11:36:35 2010 * * @brief Read/Write for MSH files generated by gmsh * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* ----------------------------------------------------------------------------- Version (Legacy) 1.0 $NOD number-of-nodes node-number x-coord y-coord z-coord ... $ENDNOD $ELM number-of-elements elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list ... $ENDELM ----------------------------------------------------------------------------- Version 2.1 $MeshFormat version-number file-type data-size $EndMeshFormat $Nodes number-of-nodes node-number x-coord y-coord z-coord ... $EndNodes $Elements number-of-elements elm-number elm-type number-of-tags < tag > ... node-number-list ... $EndElements $PhysicalNames number-of-names physical-dimension physical-number "physical-name" ... $EndPhysicalNames $NodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... node-number value ... ... $EndNodeData $ElementData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number value ... ... $EndElementData $ElementNodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number number-of-nodes-per-element value ... ... $ElementEndNodeData ----------------------------------------------------------------------------- elem-type 1: 2-node line. 2: 3-node triangle. 3: 4-node quadrangle. 4: 4-node tetrahedron. 5: 8-node hexahedron. 6: 6-node prism. 7: 5-node pyramid. 8: 3-node second order line 9: 6-node second order triangle 10: 9-node second order quadrangle 11: 10-node second order tetrahedron 12: 27-node second order hexahedron 13: 18-node second order prism 14: 14-node second order pyramid 15: 1-node point. 16: 8-node second order quadrangle 17: 20-node second order hexahedron 18: 15-node second order prism 19: 13-node second order pyramid 20: 9-node third order incomplete triangle 21: 10-node third order triangle 22: 12-node fourth order incomplete triangle 23: 15-node fourth order triangle 24: 15-node fifth order incomplete triangle 25: 21-node fifth order complete triangle 26: 4-node third order edge 27: 5-node fourth order edge 28: 6-node fifth order edge 29: 20-node third order tetrahedron 30: 35-node fourth order tetrahedron 31: 56-node fifth order tetrahedron -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "mesh_io_msh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Constant arrays initialisation */ /* -------------------------------------------------------------------------- */ UInt MeshIOMSH::_read_order[_max_element_type][MAX_NUMBER_OF_NODE_PER_ELEMENT] = { { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // _not_defined { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // _line1 { 0, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, // _line2 { 0, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, // _triangle_3 { 0, 1, 2, 3, 4, 5, 6, 0, 0, 0 }, // _triangle_6 { 0, 1, 2, 3, 4, 0, 0, 0, 0, 0 }, // _tetrahedron_4 { 0, 1, 2, 3, 4, 5, 6, 7, 9, 8 }, // _tetrahedron_10 { 0, 1, 2, 3, 0, 0, 0, 0, 0, 0 }, // _quadrangle_4 - { 0, 1, 2, 3, 4, 5, 6, 7, 0, 0 } // _hexahedron_8 + { 0, 1, 2, 3, 4, 5, 6, 7, 0, 0 }, // _quadrangle_6 + { 0, 1, 2, 3, 4, 5, 6, 7, 0, 0 } // _hexahedron_8 }; -UInt MeshIOMSH::_msh_nodes_per_elem[16] = +UInt MeshIOMSH::_msh_nodes_per_elem[17] = { 0, // element types began at 1 2, 3, 4, 4, 8, 6, 5, // 1st order 3, 6, 9, 10, 27, 18, 14, // 2st order - 1 + 1, 8 }; -ElementType MeshIOMSH::_msh_to_akantu_element_types[16] = +ElementType MeshIOMSH::_msh_to_akantu_element_types[17] = { _not_defined, // element types began at 1 _segment_2, _triangle_3, _quadrangle_4, _tetrahedron_4, // 1st order _hexahedron_8, _not_defined, _not_defined, _segment_3, _triangle_6, _not_defined, _tetrahedron_10, // 2nd order _not_defined, _not_defined, _not_defined, - _not_defined + _not_defined, _quadrangle_8 }; MeshIOMSH::MSHElementType MeshIOMSH::_akantu_to_msh_element_types[_max_element_type] = { _msh_not_defined, // _not_defined _msh_segment_2, // _segment_2 _msh_segment_3, // _segment_3 _msh_triangle_3, // _triangle_3 _msh_triangle_6, // _triangle_6 _msh_tetrahedron_4, // _tetrahedron_4 _msh_tetrahedron_10,// _tetrahedron_10 - _msh_quadrangle_1, // _quadrangle_4 + _msh_quadrangle_4, // _quadrangle_4 + _msh_quadrangle_8, // _quadrangle_4 _msh_hexahedron_8 // _hexahedron_8 }; /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIOMSH::MeshIOMSH() { canReadSurface = false; canReadExtendedData = false; } /* -------------------------------------------------------------------------- */ MeshIOMSH::~MeshIOMSH() { } /* -------------------------------------------------------------------------- */ void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt first_node_number = std::numeric_limits::max(), last_node_number = 0, file_format = 1, current_line = 0; if(!infile.good()) { AKANTU_DEBUG_ERROR("Cannot open file " << filename); } while(infile.good()) { std::getline(infile, line); current_line++; /// read the header if(line == "$MeshFormat") { std::getline(infile, line); /// the format line std::stringstream sstr(line); std::string version; sstr >> version; Int format; sstr >> format; if(format != 0) AKANTU_DEBUG_ERROR("This reader can only read ASCII files."); std::getline(infile, line); /// the end of block line current_line += 2; file_format = 2; } /// read all nodes if(line == "$Nodes" || line == "$NOD") { UInt nb_nodes; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_nodes; current_line++; Vector & nodes = const_cast &>(mesh.getNodes()); nodes.resize(nb_nodes); mesh.nb_global_nodes = nb_nodes; UInt index; Real coord[3]; UInt spatial_dimension = nodes.getNbComponent(); /// for each node, read the coordinates for(UInt i = 0; i < nb_nodes; ++i) { UInt offset = i * spatial_dimension; std::getline(infile, line); std::stringstream sstr_node(line); sstr_node >> index >> coord[0] >> coord[1] >> coord[2]; current_line++; first_node_number = std::min(first_node_number,index); last_node_number = std::max(last_node_number, index); /// read the coordinates for(UInt j = 0; j < spatial_dimension; ++j) nodes.values[offset + j] = coord[j]; } std::getline(infile, line); /// the end of block line } /// read all elements if(line == "$Elements" || line == "$ELM") { UInt nb_elements; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_elements; current_line++; Int index; UInt msh_type; ElementType akantu_type, akantu_type_old = _not_defined; Vector *connectivity = NULL; UInt node_per_element = 0; for(UInt i = 0; i < nb_elements; ++i) { std::getline(infile, line); std::stringstream sstr_elem(line); current_line++; sstr_elem >> index; sstr_elem >> msh_type; /// get the connectivity vector depending on the element type akantu_type = _msh_to_akantu_element_types[msh_type]; if(akantu_type == _not_defined) continue; if(akantu_type != akantu_type_old) { connectivity = mesh.getConnectivityPointer(akantu_type); connectivity->resize(0); node_per_element = connectivity->getNbComponent(); akantu_type_old = akantu_type; } /// read tags informations if(file_format == 2) { UInt nb_tags; sstr_elem >> nb_tags; for(UInt j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; std::stringstream sstr_tag_name; sstr_tag_name << "tag_" << j; Vector * data = mesh.getUIntDataPointer(akantu_type, sstr_tag_name.str()); data->push_back(tag); } } else if (file_format == 1) { Int tag; sstr_elem >> tag; //reg-phys std::string tag_name = "tag_0"; Vector * data = mesh.getUIntDataPointer(akantu_type, tag_name); data->push_back(tag); sstr_elem >> tag; //reg-elem tag_name = "tag_1"; data = mesh.getUIntDataPointer(akantu_type, tag_name); data->push_back(tag); sstr_elem >> tag; //number-of-nodes } UInt local_connect[node_per_element]; for(UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= last_node_number, "Node number not in range : line " << current_line); node_index -= first_node_number; local_connect[_read_order[akantu_type][j]] = node_index; } connectivity->push_back(local_connect); } std::getline(infile, line); /// the end of block line } if((line[0] == '$') && (line.find("End") == std::string::npos)) { AKANTU_DEBUG_WARNING("Unsuported block_kind " << line << " at line " << current_line); } } infile.close(); } /* -------------------------------------------------------------------------- */ void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) { std::ofstream outfile; const Vector & nodes = mesh.getNodes(); outfile.open(filename.c_str()); outfile << "$MeshFormat" << std::endl; outfile << "2.1 0 8" << std::endl;; outfile << "$EndMeshFormat" << std::endl;; outfile << std::setprecision(std::numeric_limits::digits10); outfile << "$Nodes" << std::endl;; outfile << nodes.getSize() << std::endl;; outfile << std::uppercase; for(UInt i = 0; i < nodes.getSize(); ++i) { Int offset = i * nodes.getNbComponent(); outfile << i+1; for(UInt j = 0; j < nodes.getNbComponent(); ++j) { outfile << " " << nodes.values[offset + j]; } if(nodes.getNbComponent() == 2) outfile << " " << 0.; outfile << std::endl;; } outfile << std::nouppercase; outfile << "$EndNodes" << std::endl;; outfile << "$Elements" << std::endl;; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; Int nb_elements = 0; for(it = type_list.begin(); it != type_list.end(); ++it) { const Vector & connectivity = mesh.getConnectivity(*it); nb_elements += connectivity.getSize(); } outfile << nb_elements << std::endl; UInt element_idx = 1; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; const Vector & connectivity = mesh.getConnectivity(type); for(UInt i = 0; i < connectivity.getSize(); ++i) { UInt offset = i * connectivity.getNbComponent(); outfile << element_idx << " " << _akantu_to_msh_element_types[type] << " 3 1 1 0"; /// \todo write the real data in the file for(UInt j = 0; j < connectivity.getNbComponent(); ++j) { outfile << " " << connectivity.values[offset + j] + 1; } outfile << std::endl; element_idx++; } } outfile << "$EndElements" << std::endl;; outfile.close(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/mesh_utils/mesh_io/mesh_io_msh.hh b/mesh_utils/mesh_io/mesh_io_msh.hh index 40e350ae9..8de3606e4 100644 --- a/mesh_utils/mesh_io/mesh_io_msh.hh +++ b/mesh_utils/mesh_io/mesh_io_msh.hh @@ -1,111 +1,112 @@ /** * @file mesh_io_msh.hh * @author Nicolas Richart * @date Fri Jun 18 11:30:59 2010 * * @brief Read/Write for MSH files * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_IO_MSH_HH__ #define __AKANTU_MESH_IO_MSH_HH__ /* -------------------------------------------------------------------------- */ #include "mesh_io.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class MeshIOMSH : public MeshIO { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshIOMSH(); virtual ~MeshIOMSH(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read a mesh from the file virtual void read(const std::string & filename, Mesh & mesh); /// write a mesh to a file virtual void write(const std::string & filename, const Mesh & mesh); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// MSH element types enum MSHElementType { - _msh_not_defined = 0, - _msh_segment_2 = 1, // 2-node line. - _msh_triangle_3 = 2, // 3-node triangle. - _msh_quadrangle_1 = 3, // 4-node quadrangle. - _msh_tetrahedron_4 = 4, // 4-node tetrahedron. - _msh_hexahedron_8 = 5, // 8-node hexahedron. - _msh_prism_1 = 6, // 6-node prism. - _msh_pyramid_1 = 7, // 5-node pyramid. - _msh_segment_3 = 8, // 3-node second order line - _msh_triangle_6 = 9, // 6-node second order triangle - _msh_quadrangle_2 = 10, // 9-node second order quadrangle - _msh_tetrahedron_10 = 11, // 10-node second order tetrahedron - _msh_hexaedron_2 = 12, // 27-node second order hexahedron - _msh_prism_2 = 13, // 18-node second order prism - _msh_pyramid_2 = 14, // 14-node second order pyramid - _msh_point = 15 // 1-node point. + _msh_not_defined = 0, + _msh_segment_2 = 1, // 2-node line. + _msh_triangle_3 = 2, // 3-node triangle. + _msh_quadrangle_4 = 3, // 4-node quadrangle. + _msh_tetrahedron_4 = 4, // 4-node tetrahedron. + _msh_hexahedron_8 = 5, // 8-node hexahedron. + _msh_prism_1 = 6, // 6-node prism. + _msh_pyramid_1 = 7, // 5-node pyramid. + _msh_segment_3 = 8, // 3-node second order line + _msh_triangle_6 = 9, // 6-node second order triangle + _msh_quadrangle_9 = 10, // 9-node second order quadrangle + _msh_tetrahedron_10 = 11, // 10-node second order tetrahedron + _msh_hexaedron_2 = 12, // 27-node second order hexahedron + _msh_prism_2 = 13, // 18-node second order prism + _msh_pyramid_2 = 14, // 14-node second order pyramid + _msh_point = 15, // 1-node point. + _msh_quadrangle_8 = 16 // 7-node second order quadrangle }; #define MAX_NUMBER_OF_NODE_PER_ELEMENT 10 // tetrahedron of second order /// order in witch element as to be read static UInt _read_order[_max_element_type][MAX_NUMBER_OF_NODE_PER_ELEMENT]; /// number of nodes per msh element - static UInt _msh_nodes_per_elem[16]; // 16 = number of recognized - // msh element types +1 (for 0) + static UInt _msh_nodes_per_elem[17]; // 17 = number of recognized + // msh element types +1 (for 0) /// correspondance between msh element types and akantu element types - static ElementType _msh_to_akantu_element_types[16]; + static ElementType _msh_to_akantu_element_types[17]; /// correspondance between akantu element types and msh element types static MSHElementType _akantu_to_msh_element_types[_max_element_type]; }; __END_AKANTU__ #endif /* __AKANTU_MESH_IO_MSH_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/beam_3d.geo b/test/test_model/test_solid_mechanics_model/beam_3d.geo index d26a560de..59483264f 100644 --- a/test/test_model/test_solid_mechanics_model/beam_3d.geo +++ b/test/test_model/test_solid_mechanics_model/beam_3d.geo @@ -1,43 +1,43 @@ -cl1 = 0.1; +cl1 = 0.3; Point(1) = {0, 0, 0, cl1}; Point(2) = {10, 0, 0, cl1}; Point(3) = {10, 1, 0, cl1}; Point(4) = {5, 1, 0, cl1}; Point(5) = {0, 1, 0, cl1}; Point(6) = {0, 0, 1, cl1}; Point(7) = {10, 0, 1, cl1}; Point(11) = {10, 1, 1, cl1}; Point(15) = {5, 1, 1, cl1}; Point(19) = {0, 1, 1, cl1}; Line(1) = {1, 2}; Line(2) = {2, 3}; Line(3) = {3, 4}; Line(4) = {4, 5}; Line(5) = {5, 1}; Line(10) = {6, 7}; Line(11) = {7, 11}; Line(12) = {11, 15}; Line(13) = {15, 19}; Line(14) = {19, 6}; Line(16) = {1, 6}; Line(17) = {2, 7}; Line(21) = {3, 11}; Line(25) = {4, 15}; Line(29) = {5, 19}; Line Loop(7) = {1, 2, 3, 4, 5}; Plane Surface(7) = {7}; Line Loop(18) = {1, 17, -10, -16}; Ruled Surface(18) = {18}; Line Loop(22) = {2, 21, -11, -17}; Ruled Surface(22) = {22}; Line Loop(26) = {3, 25, -12, -21}; Ruled Surface(26) = {26}; Line Loop(30) = {4, 29, -13, -25}; Ruled Surface(30) = {30}; Line Loop(34) = {5, 16, -14, -29}; Ruled Surface(34) = {34}; Line Loop(35) = {10, 11, 12, 13, 14}; Plane Surface(35) = {35}; Surface Loop(37) = {7, 18, 22, 26, 30, 34, 35}; Volume(37) = {37}; Physical Volume(1) = {37}; diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/patch_tests/CMakeLists.txt index 980d0631c..ddc96f7da 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/patch_tests/CMakeLists.txt @@ -1,5 +1,17 @@ +set(LIST_TYPES + segment_2 + segment_3 + triangle_3 + triangle_6 + quadrangle_4 + quadrangle_8 + tetrahedron_4 + tetrahedron_10 + hexahedron_8 + ) + add_subdirectory(explicit) if(AKANTU_MUMPS_ON) add_subdirectory(implicit) endif() \ No newline at end of file diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/CMakeLists.txt index af90dac12..39e852d5c 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/CMakeLists.txt @@ -1,48 +1,37 @@ #=============================================================================== # @file CMakeLists.txt # @author David Kammer # @date Wed Feb 16 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== -set(LIST_TYPES - segment_2 - segment_3 - triangle_3 - triangle_6 - quadrangle_4 - tetrahedron_4 - tetrahedron_10 - hexahedron_8 - ) - foreach(_type ${LIST_TYPES}) register_test(patch_test_explicit_${_type} patch_test_explicit.cc) set_target_properties(patch_test_explicit_${_type} PROPERTIES COMPILE_DEFINITIONS TYPE=_${_type}) file(COPY ../data/${_type}.msh DESTINATION .) endforeach() #=============================================================================== file(COPY ../data/material_check_stress.dat DESTINATION .) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) \ No newline at end of file diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc index f145a139c..70b29881c 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc +++ b/test/test_model/test_solid_mechanics_model/patch_tests/explicit/patch_test_explicit.cc @@ -1,370 +1,371 @@ /** * @file test_check_stress.cc * @author David Kammer * @date Wed Feb 16 13:56:42 2011 * * @brief patch test for elastic material in solid mechanics model * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "element_class.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; Real alpha [3][4] = { { 0.01, 0.02, 0.03, 0.04 }, { 0.05, 0.06, 0.07, 0.08 }, { 0.09, 0.10, 0.11, 0.12 } }; #ifdef AKANTU_USE_IOHELPER void paraviewInit(Dumper & dumper, const SolidMechanicsModel & model); void paraviewDump(Dumper & dumper); #endif /* -------------------------------------------------------------------------- */ template types::Matrix prescribed_strain() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix strain(spatial_dimension, spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { strain(i, j) = alpha[i][j + 1]; } } return strain; } template types::Matrix prescribed_stress() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix stress(spatial_dimension, spatial_dimension); //plane strain in 2d types::Matrix strain(spatial_dimension, spatial_dimension); types::Matrix pstrain; pstrain = prescribed_strain(); Real nu = 0.3; Real E = 2.1e11; Real trace = 0; /// symetric part of the strain tensor for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); Real Ep = E / (1 + nu); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { stress(i, j) = Ep * strain(i,j); if(i == j) stress(i, j) += Ep * (nu / (1 - 2*nu)) * trace; } return stress; } /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { debug::setDebugLevel(dblWarning); UInt dim = ElementClass::getSpatialDimension(); const ElementType element_type = TYPE; // UInt imposing_steps = 1000; // Real max_displacement = -0.01; UInt damping_steps = 400000; UInt damping_interval = 50; Real damping_ratio = 0.99; UInt additional_steps = 20000; UInt max_steps = damping_steps + additional_steps; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; std::stringstream filename; filename << TYPE << ".msh"; mesh_io.read(filename.str(), my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors my_model.getForce().clear(); my_model.getVelocity().clear(); my_model.getAcceleration().clear(); my_model.getDisplacement().clear(); my_model.initModel(); my_model.readMaterials("material_check_stress.dat"); my_model.initMaterials(); Real time_step = my_model.getStableTimeStep()/5.; my_model.setTimeStep(time_step); my_model.assembleMassLumped(); std::cout << "The number of time steps is: " << max_steps << " (" << time_step << "s)" << std::endl; // // boundary conditions // Vector top_nodes(0, 1); // Vector base_nodes(0, 1); const Vector & coordinates = my_mesh.getNodes(); Vector & displacement = my_model.getDisplacement(); Vector & boundary = my_model.getBoundary(); MeshUtils::buildFacets(my_mesh); MeshUtils::buildSurfaceID(my_mesh); CSR surface_nodes; MeshUtils::buildNodesPerSurface(my_mesh, surface_nodes); CSR::iterator snode = surface_nodes.begin(0); for (UInt s = 0; s < surface_nodes.getNbRows(); ++s) { CSR::iterator snode = surface_nodes.begin(s); for(; snode != surface_nodes.end(s); ++snode) { UInt n = *snode; std::cout << "Node " << n << std::endl; for (UInt i = 0; i < dim; ++i) { displacement(n, i) = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); } boundary(n, i) = true; } // if (coordinates(n, 0) < Math::tolerance) { // boundary(n, 0) = true; // displacement(n, 0) = 0.; // } // if (coordinates(n, 1) < Math::tolerance) { // boundary(n, 1) = true; // displacement(n, 1) = 0.; // } } } Vector & velocity = my_model.getVelocity(); #ifdef AKANTU_USE_IOHELPER my_model.updateResidual(); DumperParaview dumper; paraviewInit(dumper, my_model); #endif std::ofstream energy; std::stringstream energy_filename; energy_filename << "energy_" << TYPE << ".csv"; energy.open(energy_filename.str().c_str()); energy << "id,time,ekin" << std::endl; Real ekin_mean = 0.; /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ for(UInt s = 1; s <= max_steps; ++s) { if(s % 10000 == 0) std::cout << "passing step " << s << "/" << max_steps << " (" << s*time_step << "s)" <(imposing_steps)) * s; // for(UInt n = 0; n < top_nodes.getSize(); ++n) { // UInt node = top_nodes(n); // displacement(node, 1) = current_displacement; // } // } // damp velocity in order to find equilibrium if((s < damping_steps) && (s % damping_interval == 0)) { velocity *= damping_ratio; } if(s % 1000 == 0) { ekin_mean = ekin_mean / 1000.; std::cout << "Ekin mean = " << ekin_mean << std::endl; if (ekin_mean < 1e-10) break; ekin_mean = 0.; } my_model.explicitPred(); my_model.updateResidual(); my_model.updateAcceleration(); my_model.explicitCorr(); akantu::Real ekin = my_model.getKineticEnergy();; ekin_mean += ekin; if(s % 1000 == 0) energy << s << "," << s*time_step << "," << ekin << std::endl; #ifdef AKANTU_USE_IOHELPER - if(s % 10000 == 0) paraviewDump(dumper); + // if(s % 10000 == 0) paraviewDump(dumper); #endif } #ifdef AKANTU_USE_IOHELPER paraviewDump(dumper); #endif energy.close(); // UInt check_element = 0; // UInt quadrature_point = 0; UInt nb_quadrature_points = ElementClass::getNbQuadraturePoint(); Vector & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); Vector & strain_vect = const_cast &>(my_model.getMaterial(0).getStrain(element_type)); Vector::iterator stress_it = stress_vect.begin(dim, dim); Vector::iterator strain_it = strain_vect.begin(dim, dim); types::Matrix presc_stress; presc_stress = prescribed_stress(); types::Matrix presc_strain; presc_strain = prescribed_strain(); UInt nb_element = my_mesh.getNbElement(TYPE); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { types::Matrix & stress = *stress_it; types::Matrix & strain = *strain_it; for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { if(!(std::abs(strain(i, j) - presc_strain(i, j)) < 1e-9)) { std::cerr << "strain[" << i << "," << j << "] = " << strain(i, j) << " but should be = " << presc_strain(i, j) << " (-" << std::abs(strain(i, j) - presc_strain(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << strain << presc_strain << std::endl; return EXIT_FAILURE; } if(!(std::abs(stress(i, j) - presc_stress(i, j)) < 1e2)) { std::cerr << "stress[" << i << "," << j << "] = " << stress(i, j) << " but should be = " << presc_stress(i, j) << " (-" << std::abs(stress(i, j) - presc_stress(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << stress << presc_stress << std::endl; return EXIT_FAILURE; } } } ++stress_it; ++strain_it; } } for (UInt n = 0; n < nb_nodes; ++n) { for (UInt i = 0; i < dim; ++i) { Real disp = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { disp += alpha[i][j + 1] * coordinates(n, j); } if(!(std::abs(displacement(n,i) - disp) < 1e-7)) { std::cerr << "displacement(" << n << ", " << i <<")=" << displacement(n,i) << " should be equal to " << disp << "(" << displacement(n,i) - disp << ")" << std::endl; return EXIT_FAILURE; } } } // std::cout << "Strain : " << strain; // std::cout << "Stress : " << stress; // finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template UInt paraviewType(); template <> UInt paraviewType<_segment_2>() { return LINE1; }; template <> UInt paraviewType<_segment_3>() { return LINE2; }; template <> UInt paraviewType<_triangle_3>() { return TRIANGLE1; }; template <> UInt paraviewType<_triangle_6>() { return TRIANGLE2; }; template <> UInt paraviewType<_quadrangle_4>() { return QUAD1; }; +template <> UInt paraviewType<_quadrangle_8>() { return QUAD2; }; template <> UInt paraviewType<_tetrahedron_4>() { return TETRA1; }; template <> UInt paraviewType<_tetrahedron_10>() { return TETRA2; }; template <> UInt paraviewType<_hexahedron_8>() { return HEX1; }; void paraviewInit(Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "out_" << TYPE; dumper.SetMode(TEXT); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getMass().values, 1, "mass"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(Dumper & dumper) { dumper.Dump(); } #endif diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/implicit/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/patch_tests/implicit/CMakeLists.txt index 1eae03706..34334a5d3 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/implicit/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/patch_tests/implicit/CMakeLists.txt @@ -1,48 +1,37 @@ #=============================================================================== # @file CMakeLists.txt # @author David Kammer # @date Wed Feb 16 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== -set(LIST_TYPES - segment_2 - segment_3 - triangle_3 - triangle_6 - quadrangle_4 - tetrahedron_4 - tetrahedron_10 - hexahedron_8 - ) - foreach(_type ${LIST_TYPES}) register_test(patch_tests_implicit_${_type} patch_test_implicit.cc) set_target_properties(patch_tests_implicit_${_type} PROPERTIES COMPILE_DEFINITIONS TYPE=_${_type}) file(COPY ../data/${_type}.msh DESTINATION .) endforeach() #=============================================================================== file(COPY ../data/material_check_stress.dat DESTINATION .) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) \ No newline at end of file diff --git a/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc b/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc index 42ef46aaa..75c28d8f6 100644 --- a/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc +++ b/test/test_model/test_solid_mechanics_model/patch_tests/implicit/patch_test_implicit.cc @@ -1,307 +1,308 @@ /** * @file test_check_stress.cc * @author David Kammer * @date Wed Feb 16 13:56:42 2011 * * @brief patch test for elastic material in solid mechanics model * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "element_class.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; Real alpha [3][4] = { { 0.01, 0.02, 0.03, 0.04 }, { 0.05, 0.06, 0.07, 0.08 }, { 0.09, 0.10, 0.11, 0.12 } }; /* -------------------------------------------------------------------------- */ template types::Matrix prescribed_strain() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix strain(spatial_dimension, spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { strain(i, j) = alpha[i][j + 1]; } } return strain; } template types::Matrix prescribed_stress() { UInt spatial_dimension = ElementClass::getSpatialDimension(); types::Matrix stress(spatial_dimension, spatial_dimension); //plane strain in 2d types::Matrix strain(spatial_dimension, spatial_dimension); types::Matrix pstrain; pstrain = prescribed_strain(); Real nu = 0.3; Real E = 2.1e11; Real trace = 0; /// symetric part of the strain tensor for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); Real Ep = E / (1 + nu); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { stress(i, j) = Ep * strain(i,j); if(i == j) stress(i, j) += Ep * (nu / (1 - 2*nu)) * trace; } return stress; } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER void paraviewInit(Dumper & dumper, const SolidMechanicsModel & model); void paraviewDump(Dumper & dumper); #endif /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { initialize(&argc, &argv); UInt dim = ElementClass::getSpatialDimension(); const ElementType element_type = TYPE; /// load mesh Mesh my_mesh(dim); MeshIOMSH mesh_io; std::stringstream filename; filename << TYPE << ".msh"; mesh_io.read(filename.str(), my_mesh); UInt nb_nodes = my_mesh.getNbNodes(); /// declaration of model SolidMechanicsModel my_model(my_mesh); /// model initialization my_model.initVectors(); // initialize the vectors my_model.getForce().clear(); my_model.getVelocity().clear(); my_model.getAcceleration().clear(); my_model.getDisplacement().clear(); my_model.initModel(); my_model.readMaterials("material_check_stress.dat"); my_model.initMaterials(); my_model.initImplicit(); const Vector & coordinates = my_mesh.getNodes(); Vector & displacement = my_model.getDisplacement(); Vector & boundary = my_model.getBoundary(); MeshUtils::buildFacets(my_mesh); MeshUtils::buildSurfaceID(my_mesh); CSR surface_nodes; MeshUtils::buildNodesPerSurface(my_mesh, surface_nodes); for (UInt s = 0; s < surface_nodes.getNbRows(); ++s) { CSR::iterator snode = surface_nodes.begin(s); for(; snode != surface_nodes.end(s); ++snode) { UInt n = *snode; std::cout << "Node " << n << std::endl; for (UInt i = 0; i < dim; ++i) { displacement(n, i) = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); } boundary(n, i) = true; } } } /* ------------------------------------------------------------------------ */ /* Main loop */ /* ------------------------------------------------------------------------ */ akantu::UInt count = 0; my_model.updateResidual(); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; paraviewInit(dumper, my_model); #endif while(!my_model.testConvergenceResidual(2e-4) && (count < 100)) { std::cout << "Iter : " << ++count << std::endl; my_model.assembleStiffnessMatrix(); my_model.solveStatic(); my_model.updateResidual(); } #ifdef AKANTU_USE_IOHELPER paraviewDump(dumper); #endif if(count > 1) { std::cerr << "The code did not converge in 1 step !" << std::endl; return EXIT_FAILURE; } /* ------------------------------------------------------------------------ */ /* Checks */ /* ------------------------------------------------------------------------ */ UInt nb_quadrature_points = ElementClass::getNbQuadraturePoint(); Vector & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); Vector & strain_vect = const_cast &>(my_model.getMaterial(0).getStrain(element_type)); Vector::iterator stress_it = stress_vect.begin(dim, dim); Vector::iterator strain_it = strain_vect.begin(dim, dim); types::Matrix presc_stress; presc_stress = prescribed_stress(); types::Matrix presc_strain; presc_strain = prescribed_strain(); UInt nb_element = my_mesh.getNbElement(TYPE); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { types::Matrix & stress = *stress_it; types::Matrix & strain = *strain_it; for (UInt i = 0; i < dim; ++i) { for (UInt j = 0; j < dim; ++j) { if(!(std::abs(strain(i, j) - presc_strain(i, j)) < 1e-15)) { std::cerr << "strain[" << i << "," << j << "] = " << strain(i, j) << " but should be = " << presc_strain(i, j) << " (-" << std::abs(strain(i, j) - presc_strain(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << "computed : " << strain << "reference : " << presc_strain << std::endl; return EXIT_FAILURE; } if(!(std::abs(stress(i, j) - presc_stress(i, j)) < 1e-3)) { std::cerr << "stress[" << i << "," << j << "] = " << stress(i, j) << " but should be = " << presc_stress(i, j) << " (-" << std::abs(stress(i, j) - presc_stress(i, j)) << ") [el : " << el<< " - q : " << q << "]" << std::endl; std::cerr << "computed : " << stress << "reference : " << presc_stress << std::endl; return EXIT_FAILURE; } } } ++stress_it; ++strain_it; } } for (UInt n = 0; n < nb_nodes; ++n) { for (UInt i = 0; i < dim; ++i) { Real disp = alpha[i][0]; for (UInt j = 0; j < dim; ++j) { disp += alpha[i][j + 1] * coordinates(n, j); } if(!(std::abs(displacement(n,i) - disp) < 1e-15)) { std::cerr << "displacement(" << n << ", " << i <<")=" << displacement(n,i) << " should be equal to " << disp << std::endl; return EXIT_FAILURE; } } } // std::cout << "Strain : " << strain; // std::cout << "Stress : " << stress; // finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER template UInt paraviewType(); template <> UInt paraviewType<_segment_2>() { return LINE1; }; template <> UInt paraviewType<_segment_3>() { return LINE2; }; template <> UInt paraviewType<_triangle_3>() { return TRIANGLE1; }; template <> UInt paraviewType<_triangle_6>() { return TRIANGLE2; }; template <> UInt paraviewType<_quadrangle_4>() { return QUAD1; }; +template <> UInt paraviewType<_quadrangle_8>() { return QUAD2; }; template <> UInt paraviewType<_tetrahedron_4>() { return TETRA1; }; template <> UInt paraviewType<_tetrahedron_10>() { return TETRA2; }; template <> UInt paraviewType<_hexahedron_8>() { return HEX1; }; void paraviewInit(Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "out_" << TYPE; dumper.SetMode(TEXT); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "force"); dumper.AddNodeDataField(model.getMass().values, 1, "mass"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(Dumper & dumper) { dumper.Dump(); } #endif diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc index d7e51ed44..c4e199308 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_implicit_dynamic_2d.cc @@ -1,287 +1,287 @@ /** * @file test_solid_mechanics_model_implicit_dynamic_2d.cc * @author Nicolas Richart * @date Fri Apr 29 11:32:25 2011 * * @brief test of the dynamic implicit code * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" #include "communicator.hh" #include "mesh_partition_scotch.hh" #ifdef AKANTU_USE_SCOTCH #include "mesh_partition_scotch.hh" #endif using namespace akantu; /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" void paraviewInit(Dumper & dumper, const SolidMechanicsModel & model); void paraviewDump(Dumper & dumper); #endif /* -------------------------------------------------------------------------- */ #define bar_length 10 #define bar_height 1 #define TYPE _tetrahedron_10 UInt spatial_dimension = 3; Real time_step = 1e-4; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { debug::setDebugLevel(dblWarning); initialize(&argc, &argv); UInt max_steps = 10000000; Mesh mesh(spatial_dimension); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm->getNbProc(); Int prank = comm->whoAmI(); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ Communicator * communicator; if(prank == 0) { MeshIOMSH mesh_io; mesh_io.read("beam_3d_quad.msh", mesh); MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->reorder(); partition->partitionate(psize); communicator = Communicator::createCommunicatorDistributeMesh(mesh, partition); delete partition; } else { communicator = Communicator::createCommunicatorDistributeMesh(mesh, NULL); } SolidMechanicsModel * model = new SolidMechanicsModel(mesh); // UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); /// model initialization model->initVectors(); model->initModel(); model->readMaterials("material.dat"); model->initMaterials(); model->registerSynchronizer(*communicator); model->initImplicit(true); model->assembleMass(); // boundary conditions const Vector & position = mesh.getNodes(); Vector & boundary = model->getBoundary(); Vector & force = model->getForce(); Vector & displacment = model->getDisplacement(); //initial conditions model->getForce().clear(); model->getDisplacement().clear(); model->getVelocity().clear(); model->getAcceleration().clear(); MeshUtils::buildFacets(mesh); MeshUtils::buildSurfaceID(mesh); CSR surface_nodes; MeshUtils::buildNodesPerSurface(mesh, surface_nodes); for (UInt s = 0; s < surface_nodes.getNbRows(); ++s) { CSR::iterator snode = surface_nodes.begin(s); for(; snode != surface_nodes.end(s); ++snode) { UInt n = *snode; Real x = position(n, 0); Real y = position(n, 1); Real z = -1; if(spatial_dimension == 3) z = position(n, 2); if(Math::are_float_equal(x, 0.) && Math::are_float_equal(y, 0.)) { boundary(n,0) = true; boundary(n,1) = true; if(spatial_dimension == 3) boundary(n,2) = true; } if(Math::are_float_equal(x, bar_length) && Math::are_float_equal(y, 0.)) { boundary(n,1) = true; if(spatial_dimension == 3) boundary(n,2) = true; } if(Math::are_float_equal(x, bar_length / 2.) && Math::are_float_equal(y, bar_height)) { force(n,1) = 10000; } } } model->setTimeStep(time_step); model->updateResidual(); // model->initialAcceleration(); std::stringstream out; out << "position_" << std::scientific << time_step << ".csv"; DumperParaview dumper; paraviewInit(dumper, *model); std::ofstream pos; pos.open(out.str().c_str()); pos << "id,time,position" << std::endl; Real time = 0; UInt count = 0; - UInt print_freq = 1; + // UInt print_freq = 1; Real error; model->getMassMatrix().saveMatrix("M.mtx"); model->assembleStiffnessMatrix(); /// time loop for (UInt s = 1; s <= max_steps && time < 0.05; ++s) { model->implicitPred(); /// convergence loop do { if(count > 0) std::cout << "passing step " << s << "/" << max_steps << " " << s * time_step << "s - " << std::setw(4) << count << " : " << std::scientific << error << "\r" << std::flush; model->updateResidual(); model->solveDynamic(); model->implicitCorr(); count++; } while(!model->testConvergenceIncrement(1e-19, error) && (count < 10000)); std::cout << "passing step " << s << "/" << max_steps << " " << s * time_step << "s - " << std::setw(4) << count << " : " << std::scientific << error << std::endl; count = 0; // if(s % print_freq == 0) { // std::cout << "passing step " << s << "/" << max_steps << " " << s * time_step << "s - " << count / print_freq << std::endl; // count = 0; // } pos << s << "," << s * time_step << "," << displacment(3,1) << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10 == 0) paraviewDump(dumper); #endif time += time_step; } delete model; finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ /* Dumper vars */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ template UInt paraviewType(); -template <> UInt paraviewType<_segment_2>() { return LINE1; }; -template <> UInt paraviewType<_segment_3>() { return LINE2; }; +template <> UInt paraviewType<_segment_2>() { return LINE1; }; +template <> UInt paraviewType<_segment_3>() { return LINE2; }; template <> UInt paraviewType<_triangle_3>() { return TRIANGLE1; }; template <> UInt paraviewType<_triangle_6>() { return TRIANGLE2; }; -template <> UInt paraviewType<_quadrangle_4>() { return QUAD1; }; -template <> UInt paraviewType<_tetrahedron_4>() { return TETRA1; }; -template <> UInt paraviewType<_tetrahedron_10>() { return TETRA2; }; -template <> UInt paraviewType<_hexahedron_8>() { return HEX1; }; +template <> UInt paraviewType<_quadrangle_4>() { return QUAD1; }; +template <> UInt paraviewType<_tetrahedron_4>() { return TETRA1; }; +template <> UInt paraviewType<_tetrahedron_10>() { return TETRA2; }; +template <> UInt paraviewType<_hexahedron_8>() { return HEX1; }; /* -------------------------------------------------------------------------- */ void paraviewInit(Dumper & dumper, const SolidMechanicsModel & model) { UInt spatial_dimension = ElementClass::getSpatialDimension(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE); std::stringstream filename; filename << "dynamic_implicit_beam_" << TYPE; dumper.SetMode(TEXT); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, spatial_dimension, nb_nodes, filename.str().c_str()); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values, paraviewType(), nb_element, C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model.getAcceleration().values, spatial_dimension, "acceleration"); dumper.AddNodeDataField(model.getResidual().values, spatial_dimension, "residual"); dumper.AddNodeDataField(model.getForce().values, spatial_dimension, "applied_force"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStrain(TYPE).values, spatial_dimension*spatial_dimension, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("applied_force", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); } /* -------------------------------------------------------------------------- */ void paraviewDump(Dumper & dumper) { dumper.Dump(); } /* -------------------------------------------------------------------------- */ #endif