diff --git a/src/fe_engine/element_class.hh b/src/fe_engine/element_class.hh index fce855a3f..5eda6c5ca 100644 --- a/src/fe_engine/element_class.hh +++ b/src/fe_engine/element_class.hh @@ -1,399 +1,399 @@ /** * @file element_class.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 21 2016 * * @brief Declaration of the ElementClass main class and the * Integration and Interpolation elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_HH__ #define __AKANTU_ELEMENT_CLASS_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /// default element class structure template struct ElementClassProperty { static const GeometricalType geometrical_type = _gt_not_defined; static const InterpolationType interpolation_type = _itp_not_defined; static const ElementKind element_kind = _ek_regular; static const UInt spatial_dimension = 0; static const GaussIntegrationType gauss_integration_type = _git_not_defined; static const UInt polynomial_degree = 0; }; /// Macro to generate the element class structures for different element types #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type = geom_type; \ static const InterpolationType interpolation_type = interp_type; \ static const ElementKind element_kind = elem_kind; \ static const UInt spatial_dimension = sp; \ static const GaussIntegrationType gauss_integration_type = \ gauss_int_type; \ static const UInt polynomial_degree = min_int_order; \ } /* -------------------------------------------------------------------------- */ /* Geometry */ /* -------------------------------------------------------------------------- */ /// Default GeometricalShape structure template struct GeometricalShape { static const GeometricalShapeType shape = _gst_point; }; /// Templated GeometricalShape with function contains template struct GeometricalShapeContains { /// Check if the point (vector in 2 and 3D) at natural coordinate coor template static inline bool contains(const vector_type & coord); }; /// Macro to generate the GeometricalShape structures for different geometrical /// types #define AKANTU_DEFINE_SHAPE(geom_type, geom_shape) \ template <> struct GeometricalShape { \ static const GeometricalShapeType shape = geom_shape; \ } /* -------------------------------------------------------------------------- */ /// Templated GeometricalElement with function getInradius template ::shape> class GeometricalElement { public: /// compute the in-radius static inline Real getInradius(__attribute__((unused)) const Matrix & coord) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// true if the natural coordinates are in the element template static inline bool contains(const vector_type & coord); public: static AKANTU_GET_MACRO_NOT_CONST(SpatialDimension, spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbNodesPerElement, nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbFacetTypes, nb_facet_types, UInt); static inline UInt getNbFacetsPerElement(UInt t); static inline UInt getNbFacetsPerElement(); static inline const MatrixProxy getFacetLocalConnectivityPerElement(UInt t = 0); protected: /// Number of nodes per element static UInt nb_nodes_per_element; /// spatial dimension of the element static UInt spatial_dimension; /// number of different facet types static UInt nb_facet_types; /// number of facets for element static UInt nb_facets[]; /// storage of the facet local connectivity static UInt facet_connectivity_vect[]; /// local connectivity of facets static UInt * facet_connectivity[]; private: /// Type of the facet elements static UInt nb_nodes_per_facet[]; }; /* -------------------------------------------------------------------------- */ /* Interpolation */ /* -------------------------------------------------------------------------- */ /// default InterpolationPorperty structure template struct InterpolationPorperty { static const InterpolationKind kind = _itk_not_defined; static const UInt nb_nodes_per_element = 0; static const UInt natural_space_dimension = 0; }; /// Macro to generate the InterpolationPorperty structures for different /// interpolation types #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) \ template <> struct InterpolationPorperty { \ static const InterpolationKind kind = itp_kind; \ static const UInt nb_nodes_per_element = nb_nodes; \ static const UInt natural_space_dimension = ndim; \ } #include "interpolation_element_tmpl.hh" /* -------------------------------------------------------------------------- */ /// Generic (templated by the enum InterpolationType which specifies the order /// and the dimension of the interpolation) class handling the elemental /// interpolation template ::kind> class InterpolationElement { public: typedef InterpolationPorperty interpolation_property; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix & natural_coord, Matrix & N); /// compute the shape values for a given point in natural coordinates template static inline void computeShapes(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) vector_type & N) { AKANTU_DEBUG_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix & natural_coord, Tensor3 & dnds); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template static inline void computeDNDS(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) matrix_type & dnds) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute jacobian (or integration variable change factor) for a given point /// in the case of spatial_dimension != natural_space_dimension static inline void computeSpecialJacobian(__attribute__((unused)) const Matrix & J, __attribute__((unused)) Real & jacobians) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate a field given (arbitrary) natural coordinates static inline void interpolateOnNaturalCoordinates(const Vector & natural_coords, const Matrix & nodal_values, Vector & interpolated); /// interpolate a field given the shape functions on the interpolation point static inline void interpolate(const Matrix & nodal_values, const Vector & shapes, Vector & interpolated); /// interpolate a field given the shape functions on the interpolations points static inline void interpolate(const Matrix & nodal_values, const Matrix & shapes, Matrix & interpolated); /// compute the gradient of a given field on the given natural coordinates static inline void gradientOnNaturalCoordinates(const Vector & natural_coords, const Matrix & f, Matrix & gradient); public: static AKANTU_GET_MACRO_NOT_CONST( ShapeSize, InterpolationPorperty::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeDerivativesSize, (InterpolationPorperty::nb_nodes_per_element * InterpolationPorperty::natural_space_dimension), UInt); static AKANTU_GET_MACRO_NOT_CONST( NaturalSpaceDimension, InterpolationPorperty::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbNodesPerInterpolationElement, InterpolationPorperty::nb_nodes_per_element, UInt); }; /* -------------------------------------------------------------------------- */ /* Integration */ /* -------------------------------------------------------------------------- */ template struct GaussIntegrationTypeData { /// quadrature points in natural coordinates static Real quad_positions[]; /// weights for the Gauss integration static Real quad_weights[]; }; template ::polynomial_degree> class GaussIntegrationElement { public: static UInt getNbQuadraturePoints(); static const Matrix getQuadraturePoints(); static const Vector getWeights(); }; /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ template ::element_kind> class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: typedef GeometricalElement< ElementClassProperty::geometrical_type> geometrical_element; typedef InterpolationElement::interpolation_type> interpolation_element; typedef ElementClassProperty element_property; typedef typename interpolation_element::interpolation_property interpolation_property; public: /** * compute @f$ J = \frac{\partial x_j}{\partial s_i} @f$ the variation of real * coordinates along with variation of natural coordinates on a given point in * natural coordinates */ static inline void computeJMat(const Matrix & dnds, const Matrix & node_coords, Matrix & J); /** * compute the Jacobian matrix by computing the variation of real coordinates * along with variation of natural coordinates on a given set of points in * natural coordinates */ static inline void computeJMat(const Tensor3 & dnds, const Matrix & node_coords, Tensor3 & J); /// compute the jacobians of a serie of natural coordinates static inline void computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians); /// compute jacobian (or integration variable change factor) for a set of /// points static inline void computeJacobian(const Tensor3 & J, Vector & jacobians); /// compute jacobian (or integration variable change factor) for a given point static inline void computeJacobian(const Matrix & J, Real & jacobians); /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3 & J, const Tensor3 & dnds, Tensor3 & shape_deriv); /// compute shape derivatives (input is dxds) for a given point static inline void computeShapeDerivatives(const Matrix & J, const Matrix & dnds, Matrix & shape_deriv); /// compute the normal of a surface defined by the function f static inline void computeNormalsOnNaturalCoordinates(const Matrix & coord, Matrix & f, Matrix & normals); /// get natural coordinates from real coordinates static inline void inverseMap(const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, - Real tolerance = 1e-8); + Real tolerance = 1e-10); /// get natural coordinates from real coordinates static inline void inverseMap(const Matrix & real_coords, const Matrix & node_coords, Matrix & natural_coords, - Real tolerance = 1e-8); + Real tolerance = 1e-10); public: static AKANTU_GET_MACRO_NOT_CONST(Kind, element_kind, ElementKind); static AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty::spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, p1_type, const ElementType &); static const ElementType & getFacetType(UInt t = 0) { return facet_type[t]; } static ElementType * getFacetTypeInternal() { return facet_type; } protected: /// Type of the facet elements static ElementType facet_type[]; /// type of element P1 associated static ElementType p1_type; }; /* -------------------------------------------------------------------------- */ } // akantu #include "element_class_tmpl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { #include "element_class_point_1_inline_impl.cc" #include "element_class_segment_2_inline_impl.cc" #include "element_class_segment_3_inline_impl.cc" #include "element_class_triangle_3_inline_impl.cc" #include "element_class_triangle_6_inline_impl.cc" #include "element_class_tetrahedron_4_inline_impl.cc" #include "element_class_tetrahedron_10_inline_impl.cc" #include "element_class_quadrangle_4_inline_impl.cc" #include "element_class_quadrangle_8_inline_impl.cc" #include "element_class_hexahedron_8_inline_impl.cc" #include "element_class_hexahedron_20_inline_impl.cc" #include "element_class_pentahedron_6_inline_impl.cc" #include "element_class_pentahedron_15_inline_impl.cc" } // akantu /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) #include "cohesive_element.hh" #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "element_class_structural.hh" #endif #if defined(AKANTU_IGFEM) #include "element_class_igfem.hh" #endif #endif /* __AKANTU_ELEMENT_CLASS_HH__ */ diff --git a/src/fe_engine/element_class_tmpl.hh b/src/fe_engine/element_class_tmpl.hh index e6bdc6db1..51e25a110 100644 --- a/src/fe_engine/element_class_tmpl.hh +++ b/src/fe_engine/element_class_tmpl.hh @@ -1,510 +1,513 @@ /** * @file element_class_tmpl.hh * * @author Aurelia Isabel Cuba Ramos * @author Thomas Menouillard * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Thu Oct 22 2015 * * @brief Implementation of the inline templated function of the element class * descriptions * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "gauss_integration_tmpl.hh" +#include "element_class.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_ELEMENT_CLASS_TMPL_HH__ +#define __AKANTU_ELEMENT_CLASS_TMPL_HH__ namespace akantu { - /* -------------------------------------------------------------------------- */ /* GeometricalElement */ /* -------------------------------------------------------------------------- */ template inline const MatrixProxy GeometricalElement::getFacetLocalConnectivityPerElement(UInt t) { return MatrixProxy(facet_connectivity[t], nb_facets[t], nb_nodes_per_facet[t]); } /* -------------------------------------------------------------------------- */ template inline UInt GeometricalElement::getNbFacetsPerElement() { UInt total_nb_facets = 0; for (UInt n = 0; n < nb_facet_types; ++n) { total_nb_facets += nb_facets[n]; } return total_nb_facets; } /* -------------------------------------------------------------------------- */ template inline UInt GeometricalElement::getNbFacetsPerElement(UInt t) { return nb_facets[t]; } /* -------------------------------------------------------------------------- */ template template inline bool GeometricalElement::contains( const vector_type & coords) { return GeometricalShapeContains::contains(coords); } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_point>::contains(const vector_type & coords) { return (coords(0) < std::numeric_limits::epsilon()); } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_square>::contains(const vector_type & coords) { bool in = true; for (UInt i = 0; i < coords.size() && in; ++i) in &= ((coords(i) >= -(1. + std::numeric_limits::epsilon())) && (coords(i) <= (1. + std::numeric_limits::epsilon()))); return in; } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_triangle>::contains(const vector_type & coords) { bool in = true; Real sum = 0; for (UInt i = 0; (i < coords.size()) && in; ++i) { in &= ((coords(i) >= -(Math::getTolerance())) && (coords(i) <= (1. + Math::getTolerance()))); sum += coords(i); } if (in) return (in && (sum <= (1. + Math::getTolerance()))); return in; } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_prism>::contains(const vector_type & coords) { bool in = ((coords(0) >= -1.) && (coords(0) <= 1.)); // x in segment [-1, 1] // y and z in triangle in &= ((coords(1) >= 0) && (coords(1) <= 1.)); in &= ((coords(2) >= 0) && (coords(2) <= 1.)); Real sum = coords(1) + coords(2); return (in && (sum <= 1)); } /* -------------------------------------------------------------------------- */ /* InterpolationElement */ /* -------------------------------------------------------------------------- */ template inline void InterpolationElement::computeShapes( const Matrix & natural_coord, Matrix & N) { UInt nb_points = natural_coord.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector Np(N(p)); Vector ncoord_p(natural_coord(p)); computeShapes(ncoord_p, Np); } } /* -------------------------------------------------------------------------- */ template inline void InterpolationElement::computeDNDS( const Matrix & natural_coord, Tensor3 & dnds) { UInt nb_points = natural_coord.cols(); for (UInt p = 0; p < nb_points; ++p) { Matrix dnds_p(dnds(p)); Vector ncoord_p(natural_coord(p)); computeDNDS(ncoord_p, dnds_p); } } /* -------------------------------------------------------------------------- */ /** * interpolate on a point a field for which values are given on the * node of the element using the shape functions at this interpolation point * * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param shapes value of shape functions at the interpolation point * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template inline void InterpolationElement::interpolate( const Matrix & nodal_values, const Vector & shapes, Vector & interpolated) { Matrix interpm(interpolated.storage(), nodal_values.rows(), 1); Matrix shapesm( shapes.storage(), InterpolationPorperty::nb_nodes_per_element, 1); interpm.mul(nodal_values, shapesm); } /* -------------------------------------------------------------------------- */ /** * interpolate on several points a field for which values are given on the * node of the element using the shape functions at the interpolation point * * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param shapes value of shape functions at the interpolation point * @param interpolated interpolated values of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template inline void InterpolationElement::interpolate( const Matrix & nodal_values, const Matrix & shapes, Matrix & interpolated) { UInt nb_points = shapes.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector Np(shapes(p)); Vector interpolated_p(interpolated(p)); interpolate(nodal_values, Np, interpolated_p); } } /* -------------------------------------------------------------------------- */ /** * interpolate the field on a point given in natural coordinates the field which * values are given on the node of the element * * @param natural_coords natural coordinates of point where to interpolate \xi * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template inline void InterpolationElement::interpolateOnNaturalCoordinates( const Vector & natural_coords, const Matrix & nodal_values, Vector & interpolated) { Vector shapes( InterpolationPorperty::nb_nodes_per_element); computeShapes(natural_coords, shapes); interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ /// @f$ gradient_{ij} = \frac{\partial f_j}{\partial s_i} = \sum_k /// \frac{\partial N_k}{\partial s_i}f_{j n_k} @f$ template inline void InterpolationElement::gradientOnNaturalCoordinates( const Vector & natural_coords, const Matrix & f, Matrix & gradient) { Matrix dnds( InterpolationPorperty::natural_space_dimension, InterpolationPorperty::nb_nodes_per_element); computeDNDS(natural_coords, dnds); gradient.mul(f, dnds); } /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJMat(const Tensor3 & dnds, const Matrix & node_coords, Tensor3 & J) { UInt nb_points = dnds.size(2); for (UInt p = 0; p < nb_points; ++p) { Matrix J_p(J(p)); Matrix dnds_p(dnds(p)); computeJMat(dnds_p, node_coords, J_p); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJMat(const Matrix & dnds, const Matrix & node_coords, Matrix & J) { /// @f$ J = dxds = dnds * x @f$ J.mul(dnds, node_coords); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians) { UInt nb_points = natural_coords.cols(); Matrix dnds(interpolation_property::natural_space_dimension, interpolation_property::nb_nodes_per_element); Matrix J(natural_coords.rows(), node_coords.rows()); for (UInt p = 0; p < nb_points; ++p) { Vector ncoord_p(natural_coords(p)); interpolation_element::computeDNDS(ncoord_p, dnds); computeJMat(dnds, node_coords, J); computeJacobian(J, jacobians(p)); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Tensor3 & J, Vector & jacobians) { UInt nb_points = J.size(2); for (UInt p = 0; p < nb_points; ++p) { computeJacobian(J(p), jacobians(p)); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Matrix & J, Real & jacobians) { if (J.rows() == J.cols()) { jacobians = Math::det(J.storage()); } else { interpolation_element::computeSpecialJacobian(J, jacobians); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Tensor3 & J, const Tensor3 & dnds, Tensor3 & shape_deriv) { UInt nb_points = J.size(2); for (UInt p = 0; p < nb_points; ++p) { Matrix shape_deriv_p(shape_deriv(p)); computeShapeDerivatives(J(p), dnds(p), shape_deriv_p); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Matrix & J, const Matrix & dnds, Matrix & shape_deriv) { Matrix inv_J(J.rows(), J.cols()); Math::inv(J.storage(), inv_J.storage()); shape_deriv.mul(inv_J, dnds); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeNormalsOnNaturalCoordinates( const Matrix & coord, Matrix & f, Matrix & normals) { UInt dimension = normals.rows(); UInt nb_points = coord.cols(); AKANTU_DEBUG_ASSERT((dimension - 1) == interpolation_property::natural_space_dimension, "cannot extract a normal because of dimension mismatch " << dimension - 1 << " " << interpolation_property::natural_space_dimension); Matrix J(dimension, interpolation_property::natural_space_dimension); for (UInt p = 0; p < nb_points; ++p) { interpolation_element::gradientOnNaturalCoordinates(coord(p), f, J); if (dimension == 2) { Math::normal2(J.storage(), normals(p).storage()); } if (dimension == 3) { Math::normal3(J(0).storage(), J(1).storage(), normals(p).storage()); } } } /* ------------------------------------------------------------------------- */ /** * In the non linear cases we need to iterate to find the natural coordinates *@f$\xi@f$ * provided real coordinates @f$x@f$. * * We want to solve: @f$ x- \phi(\xi) = 0@f$ with @f$\phi(\xi) = \sum_I N_I(\xi) *x_I@f$ * the mapping function which uses the nodal coordinates @f$x_I@f$. * * To that end we use the Newton method and the following series: * * @f$ \frac{\partial \phi(x_k)}{\partial \xi} \left( \xi_{k+1} - \xi_k \right) *= x - \phi(x_k)@f$ * * When we consider elements embedded in a dimension higher than them (2D *triangle in a 3D space for example) * @f$ J = \frac{\partial \phi(\xi_k)}{\partial \xi}@f$ is of dimension *@f$dim_{space} \times dim_{elem}@f$ which * is not invertible in most cases. Rather we can solve the problem: * * @f$ J^T J \left( \xi_{k+1} - \xi_k \right) = J^T \left( x - \phi(\xi_k) *\right) @f$ * * So that * * @f$ d\xi = \xi_{k+1} - \xi_k = (J^T J)^{-1} J^T \left( x - \phi(\xi_k) *\right) @f$ * * So that if the series converges we have: * * @f$ 0 = J^T \left( \phi(\xi_\infty) - x \right) @f$ * * And we see that this is ill-posed only if @f$ J^T x = 0@f$ which means that *the vector provided * is normal to any tangent which means it is outside of the element itself. * * @param real_coords: the real coordinates the natural coordinates are sought *for * @param node_coords: the coordinates of the nodes forming the element * @param natural_coords: output->the sought natural coordinates * @param spatial_dimension: spatial dimension of the problem * **/ template inline void ElementClass::inverseMap( const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, Real tolerance) { UInt spatial_dimension = real_coords.size(); UInt dimension = natural_coords.size(); // matrix copy of the real_coords Matrix mreal_coords(real_coords.storage(), spatial_dimension, 1); // initial guess // Matrix natural_guess(natural_coords.storage(), dimension, 1); natural_coords.clear(); // real space coordinates provided by initial guess Matrix physical_guess(dimension, 1); // objective function f = real_coords - physical_guess Matrix f(dimension, 1); // dnds computed on the natural_guess // Matrix dnds(interpolation_element::nb_nodes_per_element, // spatial_dimension); // J Jacobian matrix computed on the natural_guess Matrix J(spatial_dimension, dimension); // G = J^t * J Matrix G(spatial_dimension, spatial_dimension); // Ginv = G^{-1} Matrix Ginv(spatial_dimension, spatial_dimension); // J = Ginv * J^t Matrix F(spatial_dimension, dimension); // dxi = \xi_{k+1} - \xi in the iterative process Matrix dxi(spatial_dimension, 1); /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation - Vector physical_guess_v(physical_guess.storage(), dimension); - interpolation_element::interpolateOnNaturalCoordinates( - natural_coords, node_coords, physical_guess_v); + auto update_f = [&f, &physical_guess, &natural_coords, &node_coords, &mreal_coords, + dimension]() { + Vector physical_guess_v(physical_guess.storage(), dimension); + interpolation_element::interpolateOnNaturalCoordinates( + natural_coords, node_coords, physical_guess_v); + + // compute initial objective function value f = real_coords - physical_guess + f = mreal_coords; + f -= physical_guess; - // compute initial objective function value f = real_coords - physical_guess - f = mreal_coords; - f -= physical_guess; + // compute initial error + auto error = f.norm(); + return error; + }; - // compute initial error - Real inverse_map_error = f.norm(); + auto inverse_map_error = update_f(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ while (tolerance < inverse_map_error) { // compute J^t interpolation_element::gradientOnNaturalCoordinates(natural_coords, node_coords, J); // compute G G.mul(J, J); // inverse G Ginv.inverse(G); // compute F F.mul(Ginv, J); // compute increment dxi.mul(F, f); // update our guess natural_coords += Vector(dxi(0)); - // interpolate - interpolation_element::interpolateOnNaturalCoordinates( - natural_coords, node_coords, physical_guess_v); - - // compute error - f = mreal_coords; - f -= physical_guess; - inverse_map_error = f.norm(); + inverse_map_error = update_f(); } // memcpy(natural_coords.storage(), natural_guess.storage(), sizeof(Real) * // natural_coords.size()); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::inverseMap( const Matrix & real_coords, const Matrix & node_coords, Matrix & natural_coords, Real tolerance) { UInt nb_points = real_coords.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector X(real_coords(p)); Vector ncoord_p(natural_coords(p)); inverseMap(X, node_coords, ncoord_p, tolerance); } } -} // akantu +} // namespace akantu + +#endif /* __AKANTU_ELEMENT_CLASS_TMPL_HH__ */ diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh index dc8a0c383..fcca5d08e 100644 --- a/src/fe_engine/fe_engine_template.hh +++ b/src/fe_engine/fe_engine_template.hh @@ -1,385 +1,394 @@ /** * @file fe_engine_template.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Dec 08 2015 * * @brief templated class that calls integration and shape objects * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ +/* -------------------------------------------------------------------------- */ #ifndef __AKANTU_FE_ENGINE_TEMPLATE_HH__ #define __AKANTU_FE_ENGINE_TEMPLATE_HH__ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" #include "integrator.hh" #include "shape_functions.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { class DOFManager; namespace fe_engine { namespace details { template struct AssembleLumpedTemplateHelper; template struct AssembleFieldMatrixHelper; } // namespace details } // namespace fe_engine template struct AssembleFieldMatrixStructHelper; struct DefaultIntegrationOrderFunctor { template static inline constexpr int getOrder() { return ElementClassProperty::polynomial_degree; } }; /* -------------------------------------------------------------------------- */ template