diff --git a/common/aka_common.hh b/common/aka_common.hh index ca61fd66c..248a2c628 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,337 +1,340 @@ /** * @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 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) \ 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; \ } \ } #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) \ + (_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 + _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 _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__ #endif /* __AKANTU_COMMON_HH__ */ diff --git a/fem/element_class.cc b/fem/element_class.cc index eb84e66b2..ecee71e3f 100644 --- a/fem/element_class.cc +++ b/fem/element_class.cc @@ -1,162 +1,188 @@ /** * @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>::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<_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, 2, 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]}; +/* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/fem/element_class_inline_impl.cc b/fem/element_class_inline_impl.cc index 67c1d9810..df1bb55b0 100644 --- a/fem/element_class_inline_impl.cc +++ b/fem/element_class_inline_impl.cc @@ -1,243 +1,244 @@ /** * @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 computeJacobian(dxds, nb_quadrature_points, dimension, jacobians); // if dimension == spatial_dimension compute shape derivatives if (dimension == spatial_dimension) { computeShapeDerivatives(dxds, dnds, nb_quadrature_points, dimension, dshape); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes) { Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapes(cpoint, shapes); shapes += nb_nodes_per_element; cpoint += spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(const Real * natural_coords, const UInt nb_points, Real * dnds) { Real * cpoint = const_cast(natural_coords); Real * cdnds = dnds; for (UInt p = 0; p < nb_points; ++p) { computeDNDS(cpoint, cdnds); cpoint += spatial_dimension; cdnds += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const UInt nb_points, const Real * node_coords, const UInt dimension, Real * dxds) { Real * cdnds = const_cast(dnds); Real * cdxds = dxds; for (UInt p = 0; p < nb_points; ++p) { computeDXDS(cdnds, node_coords, dimension, cdxds); cdnds += nb_nodes_per_element * spatial_dimension; cdxds += spatial_dimension * dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const Real * node_coords, const UInt dimension, Real * dxds) { /// @f$ J = dxds = dnds * x @f$ Math::matrix_matrix(spatial_dimension, dimension, nb_nodes_per_element, dnds, node_coords, dxds); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Real * dxds, const UInt nb_points, const UInt dimension, Real * jac) { Real * cdxds = const_cast(dxds); Real * cjac = jac; for (UInt p = 0; p < nb_points; ++p) { computeJacobian(cdxds, dimension, *cjac); AKANTU_DEBUG_ASSERT((cjac[0] > 0), "Negative jacobian computed, possible problem in the element node order."); cdxds += spatial_dimension * dimension; cjac++; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, const UInt nb_points, const UInt dimension, Real * shape_deriv) { AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"gradient in space " << dimension << " cannot be evaluated for element of dimension " << spatial_dimension); Real * cdxds = const_cast(dxds); Real * cdnds = const_cast(dnds); for (UInt p = 0; p < nb_points; ++p) { computeShapeDerivatives(cdxds, cdnds, shape_deriv); cdnds += spatial_dimension * nb_nodes_per_element; cdxds += spatial_dimension * spatial_dimension; shape_deriv += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, Real * shape_deriv) { /// @f$ dxds = J^{-1} @f$ Real inv_dxds[spatial_dimension * spatial_dimension]; if (spatial_dimension == 1) inv_dxds[0] = 1./dxds[0]; if (spatial_dimension == 2) Math::inv2(dxds, inv_dxds); if (spatial_dimension == 3) Math::inv3(dxds, inv_dxds); Math::matrixt_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, dnds, inv_dxds, shape_deriv); } /* -------------------------------------------------------------------------- */ template inline Real ElementClass::getInradius(__attribute__ ((unused)) const Real * coord) { AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); return 0; } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeNormalsOnQuadPoint(const Real * coord, const UInt dimension, Real * normals) { AKANTU_DEBUG_ASSERT((dimension - 1) == spatial_dimension, "cannot extract a normal because of dimension mismatch " << dimension << " " << spatial_dimension); Real * cpoint = const_cast(quad); Real * cnormals = normals; Real dnds[spatial_dimension*nb_nodes_per_element]; Real dxds[spatial_dimension*dimension]; for (UInt p = 0; p < nb_quadrature_points; ++p) { computeDNDS(cpoint,dnds); computeDXDS(dnds,coord,dimension,dxds); if (dimension == 2) { Math::normal2(dxds,cnormals); } if (dimension == 3){ Math::normal3(dxds,dxds+dimension,cnormals); } cpoint += spatial_dimension; cnormals += dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * shapes){ AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * dnds){ AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(__attribute__ ((unused)) const Real * dxds, __attribute__ ((unused)) const UInt dimension, __attribute__ ((unused)) Real & jac){ AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } #include "element_classes/element_class_segment_2.cc" #include "element_classes/element_class_segment_3.cc" #include "element_classes/element_class_triangle_3.cc" #include "element_classes/element_class_triangle_6.cc" #include "element_classes/element_class_tetrahedron_4.cc" #include "element_classes/element_class_tetrahedron_10.cc" #include "element_classes/element_class_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 new file mode 100644 index 000000000..685999506 --- /dev/null +++ b/fem/element_classes/element_class_hexahedron_8.cc @@ -0,0 +1,211 @@ +/** + * @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; + } 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 015ff2c28..7b3aa683f 100644 --- a/fem/element_classes/element_class_quadrangle_4.cc +++ b/fem/element_classes/element_class_quadrangle_4.cc @@ -1,144 +1,145 @@ /** * @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} + * @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; } 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_triangle_6.cc b/fem/element_classes/element_class_triangle_6.cc index 4a91490a9..ea17a1ebc 100644 --- a/fem/element_classes/element_class_triangle_6.cc +++ b/fem/element_classes/element_class_triangle_6.cc @@ -1,207 +1,208 @@ /** * @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 N1}{\partial \xi} = - 1 + 4 \xi - * & \frac{\partial N1}{\partial \eta} = 0 \\ + * & \frac{\partial N2}{\partial \xi} = - 1 + 4 \xi + * & \frac{\partial N2}{\partial \eta} = 0 \\ * N3 = - \eta (1 - 2 \eta) - * & \frac{\partial N1}{\partial \xi} = 0 - * & \frac{\partial N1}{\partial \eta} = - 1 + 4 \eta \\ + * & \frac{\partial N3}{\partial \xi} = 0 + * & \frac{\partial N3}{\partial \eta} = - 1 + 4 \eta \\ * N4 = 4 \xi (1 - \xi - \eta) - * & \frac{\partial N1}{\partial \xi} = 4 (1 - 2 \xi - \eta) - * & \frac{\partial N1}{\partial \eta} = - 4 \eta \\ + * & \frac{\partial N4}{\partial \xi} = 4 (1 - 2 \xi - \eta) + * & \frac{\partial N4}{\partial \eta} = - 4 \eta \\ * N5 = 4 \xi \eta - * & \frac{\partial N1}{\partial \xi} = 4 \xi - * & \frac{\partial N1}{\partial \eta} = 4 \eta \\ + * & \frac{\partial N5}{\partial \xi} = 4 \xi + * & \frac{\partial N5}{\partial \eta} = 4 \eta \\ * N6 = 4 \eta (1 - \xi - \eta) - * & \frac{\partial N1}{\partial \xi} = - 4 \xi - * & \frac{\partial N1}{\partial \eta} = 4 (1 - \xi - 2 \eta) + * & \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, Real * shapes){ /// Natural coordinates Real c0 = 1 - natural_coords[0] - natural_coords[1]; /// @f$ c0 = 1 - \xi - \eta @f$ Real c1 = natural_coords[0]; /// @f$ c1 = \xi @f$ Real c2 = natural_coords[1]; /// @f$ c2 = \eta @f$ shapes[0] = c0 * (2 * c0 - 1.); shapes[1] = c1 * (2 * c1 - 1.); shapes[2] = c2 * (2 * c2 - 1.); shapes[3] = 4 * c0 * c1; shapes[4] = 4 * c1 * c2; shapes[5] = 4 * c2 * c0; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_6>::computeDNDS(const Real * natural_coords, Real * dnds){ /** * @f[ * dnds = \left( * \begin{array}{cccccc} * \frac{\partial N1}{\partial \xi} * & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} * & \frac{\partial N4}{\partial \xi} * & \frac{\partial N5}{\partial \xi} * & \frac{\partial N6}{\partial \xi} \\ * * \frac{\partial N1}{\partial \eta} * & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} * & \frac{\partial N4}{\partial \eta} * & \frac{\partial N5}{\partial \eta} * & \frac{\partial N6}{\partial \eta} * \end{array} - * \right) @f] + * \right) + * @f] */ /// Natural coordinates Real c0 = 1 - natural_coords[0] - natural_coords[1]; /// @f$ c0 = 1 - \xi - \eta @f$ Real c1 = natural_coords[0]; /// @f$ c1 = \xi @f$ Real c2 = natural_coords[1]; /// @f$ c2 = \eta @f$ dnds[0] = 1 - 4 * c0; dnds[1] = 4 * c1 - 1.; dnds[2] = 0.; dnds[3] = 4 * (c0 - c1); dnds[4] = 4 * c2; dnds[5] = - 4 * c2; dnds[6] = 1 - 4 * c0; dnds[7] = 0.; dnds[8] = 4 * c2 - 1.; dnds[9] = - 4 * c1; dnds[10] = 4 * c1; dnds[11] = 4 * (c0 - c2); } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_triangle_6>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ // if element dimension is the same as the space dimension // then jacobian factor is the determinent of dxds if (dimension == spatial_dimension){ Real weight = 1./6.; jac = Math::det2(dxds)*weight; AKANTU_DEBUG_ASSERT(jac > 0, "Negative jacobian computed, possible problem in the element node order."); } else { AKANTU_DEBUG_ERROR("to implement"); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_triangle_6>::getInradius(const Real * coord) { UInt triangles[4][3] = { {0, 3, 5}, {3, 1, 4}, {3, 4, 5}, {5, 4, 2} }; Real inradius = std::numeric_limits::max(); for (UInt t = 0; t < 4; t++) { Real ir = Math::triangle_inradius(coord + triangles[t][0] * spatial_dimension, coord + triangles[t][1] * spatial_dimension, coord + triangles[t][2] * spatial_dimension); inradius = ir < inradius ? ir : inradius; } return inradius; } diff --git a/mesh_utils/mesh_io/mesh_io_msh.cc b/mesh_utils/mesh_io/mesh_io_msh.cc index c5b2260f2..5aadcf2fa 100644 --- a/mesh_utils/mesh_io/mesh_io_msh.cc +++ b/mesh_utils/mesh_io/mesh_io_msh.cc @@ -1,416 +1,418 @@ /** * @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 }; UInt MeshIOMSH::_msh_nodes_per_elem[16] = { 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 }; ElementType MeshIOMSH::_msh_to_akantu_element_types[16] = { _not_defined, // element types began at 1 _segment_2, _triangle_3, _quadrangle_4, _tetrahedron_4, // 1st order - _not_defined, _not_defined, _not_defined, + _hexahedron_8, _not_defined, _not_defined, _segment_3, _triangle_6, _not_defined, _tetrahedron_10, // 2nd order _not_defined, _not_defined, _not_defined, _not_defined }; 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_1, // _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 = 0, 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::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 = first_node_number < index ? first_node_number : index; last_node_number = last_node_number > index ? 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; ///@todo read to get extended information on elements } } else if (file_format == 1) { Int tag; sstr_elem >> tag; //reg-phys sstr_elem >> tag; //reg-elem 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 + 1; 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"; 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 0eaf0b06c..40e350ae9 100644 --- a/mesh_utils/mesh_io/mesh_io_msh.hh +++ b/mesh_utils/mesh_io/mesh_io_msh.hh @@ -1,111 +1,111 @@ /** * @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_hexaedron_1 = 5, // 8-node hexahedron. + _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. }; #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) /// correspondance between msh element types and akantu element types static ElementType _msh_to_akantu_element_types[16]; /// 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_fem/CMakeLists.txt b/test/test_fem/CMakeLists.txt index 491d45282..498e6edca 100644 --- a/test/test_fem/CMakeLists.txt +++ b/test/test_fem/CMakeLists.txt @@ -1,130 +1,142 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # 1D #=============================================================================== add_mesh(test_fem_line_1_mesh line.geo 1 1 OUTPUT line1.msh) add_mesh(test_fem_line_2_mesh line.geo 1 2 OUTPUT line2.msh) #=============================================================================== # 1st order register_test(test_integrate_segment_2 test_integrate_segment_2.cc) add_dependencies(test_integrate_segment_2 test_fem_line_1_mesh) #=============================================================================== register_test(test_interpolate_segment_2 test_interpolate_segment_2.cc) add_dependencies(test_interpolate_segment_2 test_fem_line_1_mesh) #=============================================================================== register_test(test_gradient_segment_2 test_gradient_segment_2.cc) add_dependencies(test_gradient_segment_2 test_fem_line_1_mesh) #=============================================================================== # 2nd order register_test(test_integrate_segment_3 test_integrate_segment_3.cc) add_dependencies(test_integrate_segment_3 test_fem_line_2_mesh) #=============================================================================== register_test(test_interpolate_segment_3 test_interpolate_segment_3.cc) add_dependencies(test_interpolate_segment_3 test_fem_line_2_mesh) #=============================================================================== register_test(test_gradient_segment_3 test_gradient_segment_3.cc) add_dependencies(test_gradient_segment_3 test_fem_line_2_mesh) #=============================================================================== # 2D #=============================================================================== add_mesh(test_fem_square_1_mesh square.geo 2 1 OUTPUT square1.msh) add_mesh(test_fem_square_2_mesh square.geo 2 2 OUTPUT square2.msh) add_mesh(test_fem_circle_1_mesh circle.geo 2 1 OUTPUT circle1.msh) add_mesh(test_fem_circle_2_mesh circle.geo 2 2 OUTPUT circle2.msh) add_mesh(test_fem_square_struct_1_mesh square_structured.geo 2 1 OUTPUT square_structured1.msh) #=============================================================================== # 1st order unstructured register_test(test_integrate_triangle_3 test_integrate_triangle_3.cc) add_dependencies(test_integrate_triangle_3 test_fem_square_1_mesh test_fem_circle_1_mesh) #=============================================================================== register_test(test_interpolate_triangle_3 test_interpolate_triangle_3.cc) add_dependencies(test_interpolate_triangle_3 test_fem_square_1_mesh) #=============================================================================== register_test(test_gradient_triangle_3 test_gradient_triangle_3.cc) add_dependencies(test_gradient_triangle_3 test_fem_square_1_mesh) #=============================================================================== # 1st order structured register_test(test_integrate_quadrangle_4 test_integrate_quadrangle_4.cc) add_dependencies(test_integrate_quadrangle_4 test_fem_square_1_mesh test_fem_circle_1_mesh) #=============================================================================== register_test(test_interpolate_quadrangle_4 test_interpolate_quadrangle_4.cc) add_dependencies(test_interpolate_quadrangle_4 test_fem_square_struct_1_mesh) #=============================================================================== register_test(test_gradient_quadrangle_4 test_gradient_quadrangle_4.cc) add_dependencies(test_gradient_quadrangle_4 test_fem_square_struct_1_mesh) #=============================================================================== # 2nd order register_test(test_integrate_triangle_6 test_integrate_triangle_6.cc) add_dependencies(test_integrate_triangle_6 test_fem_square_2_mesh test_fem_circle_2_mesh) #=============================================================================== register_test(test_interpolate_triangle_6 test_interpolate_triangle_6.cc) add_dependencies(test_interpolate_triangle_6 test_fem_square_2_mesh) #=============================================================================== register_test(test_gradient_triangle_6 test_gradient_triangle_6.cc) add_dependencies(test_gradient_triangle_6 test_fem_square_2_mesh) #=============================================================================== # 3D #=============================================================================== add_mesh(test_fem_cube_1_mesh cube.geo 3 1 OUTPUT cube1.msh) add_mesh(test_fem_cube_2_mesh cube.geo 3 2 OUTPUT cube2.msh) +add_mesh(test_fem_hexa_struct_1_mesh hexa_structured.geo 3 1 OUTPUT hexa_structured1.msh) + +#=============================================================================== +# 1st order structured +register_test(test_integrate_hexahedron_8 test_integrate_hexahedron_8.cc) +add_dependencies(test_integrate_hexahedron_8 test_fem_hexa_struct_1_mesh) +#=============================================================================== +register_test(test_interpolate_hexahedron_8 test_interpolate_hexahedron_8.cc) +add_dependencies(test_interpolate_hexahedron_8 test_fem_hexa_struct_1_mesh) +#=============================================================================== +register_test(test_gradient_hexahedron_8 test_gradient_hexahedron_8.cc) +add_dependencies(test_gradient_hexahedron_8 test_fem_hexa_struct_1_mesh) #=============================================================================== # 1st order register_test(test_integrate_tetrahedron_4 test_integrate_tetrahedron_4.cc) add_dependencies(test_integrate_tetrahedron_4 test_fem_cube_1_mesh) #=============================================================================== register_test(test_interpolate_tetrahedron_4 test_interpolate_tetrahedron_4.cc) add_dependencies(test_interpolate_tetrahedron_4 test_fem_cube_1_mesh) #=============================================================================== register_test(test_gradient_tetrahedron_4 test_gradient_tetrahedron_4.cc) add_dependencies(test_gradient_tetrahedron_4 test_fem_cube_1_mesh) #=============================================================================== # 2nd order register_test(test_integrate_tetrahedron_10 test_integrate_tetrahedron_10.cc) add_dependencies(test_integrate_tetrahedron_10 test_fem_cube_2_mesh) #=============================================================================== register_test(test_interpolate_tetrahedron_10 test_interpolate_tetrahedron_10.cc) add_dependencies(test_interpolate_tetrahedron_10 test_fem_cube_2_mesh) #=============================================================================== register_test(test_gradient_tetrahedron_10 test_gradient_tetrahedron_10.cc) add_dependencies(test_gradient_tetrahedron_10 test_fem_cube_2_mesh) diff --git a/test/test_fem/hexa_structured.geo b/test/test_fem/hexa_structured.geo new file mode 100644 index 000000000..0cc869dac --- /dev/null +++ b/test/test_fem/hexa_structured.geo @@ -0,0 +1,40 @@ +h = 0.25; + +Point(1) = {0.0, 0.0, 0.0, h}; +Point(2) = {1.0, 0.0, 0.0, h}; +Point(3) = {0.0, 1.0, 0.0, h}; +Point(4) = {1.0, 1.0, 0.0, h}; +Point(5) = {0.0, 0.0, 1.0, h}; +Point(6) = {1.0, 0.0, 1.0, h}; +Point(7) = {0.0, 1.0, 1.0, h}; +Point(8) = {1.0, 1.0, 1.0, h}; + +Line(1) = {7, 8}; +Line(2) = {8, 4}; +Line(3) = {4, 3}; +Line(4) = {3, 7}; +Line(5) = {1, 5}; +Line(6) = {5, 6}; +Line(7) = {6, 2}; +Line(8) = {2, 1}; +Line(9) = {3, 1}; +Line(10) = {7, 5}; +Line(11) = {8, 6}; +Line(12) = {4, 2}; +Line Loop(13) = {1, 11, -6, -10}; +Plane Surface(14) = {13}; +Line Loop(15) = {3, 4, 1, 2}; +Plane Surface(16) = {15}; +Line Loop(17) = {6, 7, 8, 5}; +Plane Surface(18) = {17}; +Line Loop(19) = {3, 9, -8, -12}; +Plane Surface(20) = {19}; +Line Loop(21) = {4, 10, -5, -9}; +Plane Surface(22) = {21}; +Line Loop(23) = {11, 7, -12, -2}; +Plane Surface(24) = {23}; +Surface Loop(25) = {24, 14, 16, 20, 22, 18}; +Volume(26) = {25}; +Transfinite Surface "*"; +Recombine Surface "*"; +Transfinite Volume "*"; diff --git a/test/test_fem/test_gradient_XXXX.cc b/test/test_fem/test_gradient_XXXX.cc index 7d2f1316c..445677e4f 100644 --- a/test/test_fem/test_gradient_XXXX.cc +++ b/test/test_fem/test_gradient_XXXX.cc @@ -1,100 +1,104 @@ /** * @file test_gradient_XXXX.cc * @author Nicolas Richart * @date Mon Jul 19 10:55:49 2010 * * @brief test of the fem class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "fem.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "shape_lagrange.hh" #include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { ElementType type = XXXX; UInt dim = DIM; MeshIOMSH mesh_io; Mesh my_mesh(dim); mesh_io.read("FILE.msh", my_mesh); - FEMTemplate fem(my_mesh, dim, "my_fem"); + FEM *fem = new FEMTemplate(my_mesh, dim, "my_fem"); debug::setDebugLevel(dblDump); - fem.initShapeFunctions(); + fem->initShapeFunctions(); std::cout << fem << std::endl; - Vector const_val(fem.getMesh().getNbNodes(), 2, "const_val"); + StaticMemory * st_mem = StaticMemory::getStaticMemory(); + std::cout << *st_mem << std::endl; + + Vector const_val(fem->getMesh().getNbNodes(), 2, "const_val"); Vector grad_on_quad(0, 2 * dim, "grad_on_quad"); for (UInt i = 0; i < const_val.getSize(); ++i) { const_val.values[i * 2 + 0] = 1.; const_val.values[i * 2 + 1] = 2.; } - fem.gradientOnQuadraturePoints(const_val, grad_on_quad, 2, type); + fem->gradientOnQuadraturePoints(const_val, grad_on_quad, 2, type); std::ofstream my_file("out.txt"); my_file << const_val << std::endl; my_file << grad_on_quad << std::endl; // compute gradient of coordinates Vector grad_coord_on_quad(0, dim * dim, "grad_coord_on_quad"); - fem.gradientOnQuadraturePoints(my_mesh.getNodes(), grad_coord_on_quad, my_mesh.getSpatialDimension(), type); + fem->gradientOnQuadraturePoints(my_mesh.getNodes(), grad_coord_on_quad, my_mesh.getSpatialDimension(), type); my_file << my_mesh.getNodes() << std::endl; my_file << grad_coord_on_quad << std::endl; - UInt nb_quads = my_mesh.getNbElement(type) * FEM::getNbQuadraturePoints(type); - Real eps = 25 * std::numeric_limits::epsilon(); - std::cout << "Epsilon : " << eps << std::endl; - for (UInt q = 0; q < nb_quads; ++q) { - for (UInt i = 0; i < dim; ++i) { - for (UInt j = 0; j < dim; ++j) { - if(!(fabs(grad_coord_on_quad.values[q*dim*dim+ i*dim + j] - (i == j)) <= eps)) { - std::cerr << "Error on the quad point " << q << std::endl; - for (UInt oi = 0; oi < dim; ++oi) { - for (UInt oj = 0; oj < dim; ++oj) { - std::cout << fabs(grad_coord_on_quad.values[q*dim*dim + i*dim + j] - (i == j)) << " "; - } - std::cout << std::endl; - } - exit(EXIT_FAILURE); - } - } - } - } +// UInt nb_quads = my_mesh.getNbElement(type) * FEM::getNbQuadraturePoints(type); +// Real eps = 25 * std::numeric_limits::epsilon(); +// std::cout << "Epsilon : " << eps << std::endl; +// for (UInt q = 0; q < nb_quads; ++q) { +// for (UInt i = 0; i < dim; ++i) { +// for (UInt j = 0; j < dim; ++j) { +// if(!(fabs(grad_coord_on_quad.values[q*dim*dim+ i*dim + j] - (i == j)) <= eps)) { +// std::cerr << "Error on the quad point " << q << std::endl; +// for (UInt oi = 0; oi < dim; ++oi) { +// for (UInt oj = 0; oj < dim; ++oj) { +// std::cout << fabs(grad_coord_on_quad.values[q*dim*dim + i*dim + j] - (i == j)) << " "; +// } +// std::cout << std::endl; +// } +// exit(EXIT_FAILURE); +// } +// } +// } +// } + delete fem; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_fem/test_integrate_XXXX.cc b/test/test_fem/test_integrate_XXXX.cc index f10182b25..5925e2535 100644 --- a/test/test_fem/test_integrate_XXXX.cc +++ b/test/test_fem/test_integrate_XXXX.cc @@ -1,102 +1,104 @@ /** * @file test_integrate_XXXX.cc * @author Nicolas Richart * @author Guillaume Anciaux * @date Mon Jul 19 10:55:49 2010 * * @brief test of the fem class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "fem.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "shape_lagrange.hh" #include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { ElementType type = TYPE; UInt dim = DIM; MeshIOMSH mesh_io; Mesh my_mesh(dim); mesh_io.read("XXXX.msh", my_mesh); - - FEM fem(my_mesh, dim, "my_fem"); + + FEM *fem = new FEMTemplate(my_mesh, dim, "my_fem"); debug::_debug_level = dblDump; - fem.initShapeFunctions(); - + fem->initShapeFunctions(); + std::cout << *fem << std::endl; - std::cout << fem << std::endl; + StaticMemory * st_mem = StaticMemory::getStaticMemory(); + std::cout << *st_mem << std::endl; - Vector const_val(fem.getMesh().getNbNodes(), 2, "const_val"); + Vector const_val(fem->getMesh().getNbNodes(), 2, "const_val"); Vector val_on_quad(0, 2 , "val_on_quad"); for (UInt i = 0; i < const_val.getSize(); ++i) { const_val.values[i * 2 + 0] = 1.; const_val.values[i * 2 + 1] = 2.; } //interpolate function on quadrature points - fem.interpolateOnQuadraturePoints(const_val, val_on_quad, 2, type); + fem->interpolateOnQuadraturePoints(const_val, val_on_quad, 2, type); //integrate function on elements akantu::Vector int_val_on_elem(0, 2, "int_val_on_elem"); - fem.integrate(val_on_quad, int_val_on_elem, 2, type); + fem->integrate(val_on_quad, int_val_on_elem, 2, type); // get global integration value Real value[2] = {0,0}; std::ofstream my_file("out.txt"); my_file << val_on_quad << std::endl << int_val_on_elem << std::endl; - for (UInt i = 0; i < fem.getMesh().getNbElement(type); ++i) { + for (UInt i = 0; i < fem->getMesh().getNbElement(type); ++i) { value[0] += int_val_on_elem.values[2*i]; value[1] += int_val_on_elem.values[2*i+1]; } my_file << "integral is " << value[0] << " " << value[1] << std::endl; - FEM fem_boundary(my_mesh, dim-1, "my_fem_boundary"); - fem_boundary.initShapeFunctions(); + //FEM fem_boundary(my_mesh, dim-1, "my_fem_boundary"); + //fem_boundary->initShapeFunctions(); - ElementType bound_type = Mesh::getFacetElementType(type); - UInt nb_boundary_quad = FEM::getNbQuadraturePoints(bound_type); + //ElementType bound_type = Mesh::getFacetElementType(type); + //UInt nb_boundary_quad = FEM::getNbQuadraturePoints(bound_type); - Vector val_on_bquad(0, nb_boundary_quad, "val_on_quad"); - for (UInt i = 0; i < const_val.getSize(); ++i) { - const_val.values[i * 2 + 0] = 1.; - } + //Vector val_on_bquad(0, nb_boundary_quad, "val_on_quad"); + //for (UInt i = 0; i < const_val.getSize(); ++i) { + // const_val.values[i * 2 + 0] = 1.; + //} + delete fem; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_fem/test_interpolate_XXXX.cc b/test/test_fem/test_interpolate_XXXX.cc index 3297f89c2..8cfbe7193 100644 --- a/test/test_fem/test_interpolate_XXXX.cc +++ b/test/test_fem/test_interpolate_XXXX.cc @@ -1,87 +1,91 @@ /** * @file test_interpolate_XXXX.cc * @author Nicolas Richart * @date Mon Jul 19 10:55:49 2010 * * @brief test of the fem class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "fem.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "shape_lagrange.hh" #include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { ElementType type = XXXX; UInt dim = DIM; MeshIOMSH mesh_io; Mesh my_mesh(dim); mesh_io.read("XXXX.msh", my_mesh); - FEM fem(my_mesh, dim, "my_fem"); + FEM *fem = new FEMTemplate(my_mesh, dim, "my_fem"); //UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); debug::setDebugLevel(dblDump); - fem.initShapeFunctions(); + fem->initShapeFunctions(); - std::cout << fem << std::endl; + std::cout << *fem << std::endl; - Vector const_val(fem.getMesh().getNbNodes(), 2, "const_val"); + StaticMemory * st_mem = StaticMemory::getStaticMemory(); + std::cout << *st_mem << std::endl; + + Vector const_val(fem->getMesh().getNbNodes(), 2, "const_val"); Vector val_on_quad(0, 2, "val_on_quad"); for (UInt i = 0; i < const_val.getSize(); ++i) { const_val.values[i * 2 + 0] = 1.; const_val.values[i * 2 + 1] = 2.; } - fem.interpolateOnQuadraturePoints(const_val, val_on_quad, 2, type); + fem->interpolateOnQuadraturePoints(const_val, val_on_quad, 2, type); std::ofstream my_file("out.txt"); my_file << const_val << std::endl; my_file << val_on_quad << std::endl; // interpolate coordinates Vector coord_on_quad(0, my_mesh.getSpatialDimension(), "coord_on_quad"); - fem.interpolateOnQuadraturePoints(my_mesh.getNodes(), + fem->interpolateOnQuadraturePoints(my_mesh.getNodes(), coord_on_quad, my_mesh.getSpatialDimension(), type); my_file << my_mesh.getNodes() << std::endl; my_file << coord_on_quad << std::endl; + delete fem; finalize(); return EXIT_SUCCESS; }