Diffusion akantu (bugfixes/bug_fix_in_synchronization_of_mat_vec_mul_product_when_assembled_to_residual) Browse

bug_fix_in_syn…

fe_engine_template_tmpl.hh
View Options

			/**
			* @file fe_engine_template_tmpl.hh
			*
			* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
			* @author Marco Vocialta <marco.vocialta@epfl.ch>
			* @author Nicolas Richart <nicolas.richart@epfl.ch>
			*
			* @date creation: Mon Nov 05 2012
			* @date last modification: Mon Jul 07 2014
			*
			* @brief Template implementation of FEEngineTemplate
			*
			* @section LICENSE
			*
			* Copyright (©) 2014 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 <http://www.gnu.org/licenses/>.
			*
			*/

			/* -------------------------------------------------------------------------- */
			#include "aka_common.hh"

			__BEGIN_AKANTU__

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			FEEngineTemplate<I, S, kind>::FEEngineTemplate(Mesh & mesh, UInt spatial_dimension,
			ID id, MemoryID memory_id) :
			FEEngine(mesh,spatial_dimension,id,memory_id),
			integrator(mesh, id, memory_id),
			shape_functions(mesh, id, memory_id) { }

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			FEEngineTemplate<I, S, kind>::~FEEngineTemplate() { }

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct GradientOnIntegrationPointsHelper {
			template <class S>
			static void call(const S & shape_functions,
			Mesh & mesh,
			const Array<Real> & u,
			Array<Real> & nablauq,
			const UInt nb_degree_of_freedom,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define COMPUTE_GRADIENT(type) \
			if (element_dimension == ElementClass<type>::getSpatialDimension()) \
			shape_functions.template gradientOnIntegrationPoints<type>(u, \
			nablauq, \
			nb_degree_of_freedom, \
			ghost_type, \
			filter_elements);

			#define AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER(kind) \
			template<> \
			struct GradientOnIntegrationPointsHelper<kind> { \
			template <class S> \
			static void call(const S & shape_functions, \
			Mesh & mesh, \
			const Array<Real> & u, \
			Array<Real> & nablauq, \
			const UInt nb_degree_of_freedom, \
			const ElementType & type, \
			const GhostType & ghost_type, \
			const Array<UInt> & filter_elements) { \
			UInt element_dimension = mesh.getSpatialDimension(type); \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_GRADIENT, kind); \
			} \
			};
			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER, \
			AKANTU_FE_ENGINE_LIST_GRADIENT_ON_INTEGRATION_POINTS)

			#undef AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER
			#undef COMPUTE_GRADIENT

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::gradientOnIntegrationPoints(const Array<Real> &u,
			Array<Real> &nablauq,
			const UInt nb_degree_of_freedom,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) const {
			AKANTU_DEBUG_IN();

			UInt nb_element = mesh.getNbElement(type, ghost_type);
			if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
			UInt nb_points = shape_functions.getIntegrationPoints(type, ghost_type).cols();

			#ifndef AKANTU_NDEBUG

			UInt element_dimension = mesh.getSpatialDimension(type);

			AKANTU_DEBUG_ASSERT(u.getSize() == mesh.getNbNodes(),
			"The vector u(" << u.getID()
			<< ") has not the good size.");
			AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom ,
			"The vector u(" << u.getID()
			<< ") has not the good number of component.");

			AKANTU_DEBUG_ASSERT(nablauq.getNbComponent()
			== nb_degree_of_freedom * element_dimension,
			"The vector nablauq(" << nablauq.getID()
			<< ") has not the good number of component.");

			// AKANTU_DEBUG_ASSERT(nablauq.getSize() == nb_element * nb_points,
			// "The vector nablauq(" << nablauq.getID()
			// << ") has not the good size.");
			#endif

			nablauq.resize(nb_element * nb_points);

			GradientOnIntegrationPointsHelper<kind>::call(shape_functions,
			mesh,
			u,
			nablauq,
			nb_degree_of_freedom,
			type,
			ghost_type,
			filter_elements);
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::initShapeFunctions(const GhostType & ghost_type) {
			initShapeFunctions(mesh.getNodes(), ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::initShapeFunctions(const Array<Real> & nodes,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_IN();

			Mesh::type_iterator it = mesh.firstType(element_dimension, ghost_type, kind);
			Mesh::type_iterator end = mesh.lastType(element_dimension, ghost_type, kind);
			for(; it != end; ++it) {
			ElementType type = *it;
			integrator.initIntegrator(nodes, type, ghost_type);
			const Matrix<Real> & control_points =
			getIntegrationPoints(type, ghost_type);
			shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type);
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct IntegrateHelper { };

			#define INTEGRATE(type) \
			integrator.template integrate<type>(f, \
			intf, \
			nb_degree_of_freedom, \
			ghost_type, \
			filter_elements); \

			#define AKANTU_SPECIALIZE_INTEGRATE_HELPER(kind) \
			template<> \
			struct IntegrateHelper<kind> { \
			template <class I> \
			static void call(const I & integrator, \
			const Array<Real> & f, \
			Array<Real> &intf, \
			UInt nb_degree_of_freedom, \
			const ElementType & type, \
			const GhostType & ghost_type, \
			const Array<UInt> & filter_elements) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_HELPER)

			#undef AKANTU_SPECIALIZE_INTEGRATE_HELPER
			#undef INTEGRATE

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::integrate(const Array<Real> & f,
			Array<Real> &intf,
			UInt nb_degree_of_freedom,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) const{

			UInt nb_element = mesh.getNbElement(type, ghost_type);
			if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
			#ifndef AKANTU_NDEBUG

			UInt nb_quadrature_points = getNbIntegrationPoints(type);

			AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
			"The vector f(" << f.getID() << " size " << f.getSize()
			<< ") has not the good size (" << nb_element << ").");
			AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom ,
			"The vector f(" << f.getID()
			<< ") has not the good number of component.");
			AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
			"The vector intf(" << intf.getID()
			<< ") has not the good number of component.");
			#endif

			intf.resize(nb_element);

			IntegrateHelper<kind>::call(integrator,
			f,
			intf,
			nb_degree_of_freedom,
			type,
			ghost_type,
			filter_elements);

			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct IntegrateScalarHelper { };

			#define INTEGRATE(type) \
			integral = integrator.template integrate<type>(f, \
			ghost_type, filter_elements);

			#define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER(kind) \
			template<> \
			struct IntegrateScalarHelper<kind> { \
			template <class I> \
			static Real call(const I & integrator, \
			const Array<Real> & f, \
			const ElementType & type, \
			const GhostType & ghost_type, \
			const Array<UInt> & filter_elements) { \
			Real integral = 0.; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \
			return integral; \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER)

			#undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER
			#undef INTEGRATE

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			Real FEEngineTemplate<I, S, kind>::integrate(const Array<Real> & f,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) const{
			AKANTU_DEBUG_IN();

			#ifndef AKANTU_NDEBUG
			// std::stringstream sstr; sstr << ghost_type;
			// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
			// "The vector " << nablauq.getID() << " is not taged " << ghost_type);
			UInt nb_element = mesh.getNbElement(type, ghost_type);
			if(filter_elements != empty_filter) nb_element = filter_elements.getSize();

			UInt nb_quadrature_points = getNbIntegrationPoints(type, ghost_type);

			AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
			"The vector f(" << f.getID()
			<< ") has not the good size. (" << f.getSize() << "!=" << nb_quadrature_points * nb_element << ")");
			AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
			"The vector f(" << f.getID()
			<< ") has not the good number of component.");
			#endif

			Real integral = IntegrateScalarHelper<kind>::call(integrator, f, type, ghost_type, filter_elements);
			AKANTU_DEBUG_OUT();
			return integral;
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct IntegrateScalarOnOneElementHelper { };

			#define INTEGRATE(type) \
			res = integrator.template integrate<type>(f, \
			index, ghost_type);

			#define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER(kind) \
			template<> \
			struct IntegrateScalarOnOneElementHelper<kind> { \
			template <class I> \
			static Real call(const I & integrator, \
			const Vector<Real> & f, \
			const ElementType & type, \
			UInt index, \
			const GhostType & ghost_type) { \
			Real res = 0.; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \
			return res; \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER)

			#undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER
			#undef INTEGRATE

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			Real FEEngineTemplate<I, S, kind>::integrate(const Vector<Real> & f,
			const ElementType & type,
			UInt index,
			const GhostType & ghost_type) const{


			Real res = IntegrateScalarOnOneElementHelper<kind>::call(integrator, f, type, index, ghost_type);
			return res;
			}


			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct IntegrateOnIntegrationPointsHelper { };

			#define INTEGRATE(type) \
			integrator.template integrateOnQuadraturePoints<type>(f, intf, nb_degree_of_freedom, \
			ghost_type, filter_elements);

			#define AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER(kind) \
			template<> \
			struct IntegrateOnIntegrationPointsHelper<kind> { \
			template <class I> \
			static void call(const I & integrator, \
			const Array<Real> & f, \
			Array<Real> & intf, \
			UInt nb_degree_of_freedom, \
			const ElementType & type, \
			const GhostType & ghost_type, \
			const Array<UInt> & filter_elements) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER)

			#undef AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER
			#undef INTEGRATE

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::integrateOnIntegrationPoints(const Array<Real> & f,
			Array<Real> &intf,
			UInt nb_degree_of_freedom,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) const{

			UInt nb_element = mesh.getNbElement(type, ghost_type);
			if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
			UInt nb_quadrature_points = getNbIntegrationPoints(type);
			#ifndef AKANTU_NDEBUG
			// std::stringstream sstr; sstr << ghost_type;
			// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
			// "The vector " << nablauq.getID() << " is not taged " << ghost_type);


			AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
			"The vector f(" << f.getID() << " size " << f.getSize()
			<< ") has not the good size (" << nb_element << ").");
			AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom ,
			"The vector f(" << f.getID()
			<< ") has not the good number of component.");
			AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
			"The vector intf(" << intf.getID()
			<< ") has not the good number of component.");
			#endif

			intf.resize(nb_element*nb_quadrature_points);
			IntegrateOnIntegrationPointsHelper<kind>::call(integrator, f, intf, nb_degree_of_freedom, type, ghost_type, filter_elements);
			}


			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct InterpolateOnIntegrationPointsHelper {
			template <class S>
			static void call(const S & shape_functions,
			const Array<Real> &u,
			Array<Real> &uq,
			const UInt nb_degree_of_freedom,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};


			#define INTERPOLATE(type) \
			shape_functions.template interpolateOnIntegrationPoints<type>(u, \
			uq, \
			nb_degree_of_freedom, \
			ghost_type, \
			filter_elements);

			#define AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER(kind) \
			template<> \
			struct InterpolateOnIntegrationPointsHelper<kind> { \
			template <class S> \
			static void call(const S & shape_functions, \
			const Array<Real> & u, \
			Array<Real> & uq, \
			const UInt nb_degree_of_freedom, \
			const ElementType & type, \
			const GhostType & ghost_type, \
			const Array<UInt> & filter_elements) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind); \
			} \
			};
			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER, \
			AKANTU_FE_ENGINE_LIST_INTERPOLATE_ON_INTEGRATION_POINTS)

			#undef AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER
			#undef INTERPOLATE

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::interpolateOnIntegrationPoints(const Array<Real> &u,
			Array<Real> &uq,
			const UInt nb_degree_of_freedom,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) const {
			AKANTU_DEBUG_IN();

			UInt nb_points = shape_functions.getIntegrationPoints(type, ghost_type).cols();
			UInt nb_element = mesh.getNbElement(type, ghost_type);
			if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
			#ifndef AKANTU_NDEBUG

			AKANTU_DEBUG_ASSERT(u.getSize() == mesh.getNbNodes(),
			"The vector u(" << u.getID()
			<< ") has not the good size.");
			AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom ,
			"The vector u(" << u.getID()
			<< ") has not the good number of component.");
			AKANTU_DEBUG_ASSERT(uq.getNbComponent() == nb_degree_of_freedom,
			"The vector uq(" << uq.getID()
			<< ") has not the good number of component.");
			#endif

			uq.resize(nb_element * nb_points);


			InterpolateOnIntegrationPointsHelper<kind>::call(shape_functions,
			u,
			uq,
			nb_degree_of_freedom,
			type,
			ghost_type,
			filter_elements);
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::interpolateOnIntegrationPoints(const Array<Real> & u,
			ElementTypeMapArray<Real> & uq,
			const ElementTypeMapArray<UInt> * filter_elements) const {
			AKANTU_DEBUG_IN();

			const Array<UInt> * filter = NULL;

			for (ghost_type_t::iterator gt = ghost_type_t::begin();
			gt != ghost_type_t::end(); ++gt) {

			GhostType ghost_type = *gt;
			ElementTypeMapArray<Real>::type_iterator it = uq.firstType(_all_dimensions, ghost_type, kind);
			ElementTypeMapArray<Real>::type_iterator last = uq.lastType(_all_dimensions, ghost_type, kind);

			for (; it != last; ++it) {
			ElementType type = *it;

			UInt nb_quad_per_element = getNbIntegrationPoints(type, ghost_type);

			UInt nb_element = 0;

			if (filter_elements) {
			filter = &((*filter_elements)(type, ghost_type));
			nb_element = filter->getSize();
			}
			else {
			filter = &empty_filter;
			nb_element = mesh.getNbElement(type, ghost_type);
			}

			UInt nb_tot_quad = nb_quad_per_element * nb_element;

			Array<Real> & quad = uq(type, ghost_type);
			quad.resize(nb_tot_quad);

			interpolateOnIntegrationPoints(u,
			quad,
			mesh.getSpatialDimension(),
			type,
			ghost_type,
			*filter);
			}
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::computeIntegrationPointsCoordinates(ElementTypeMapArray<Real> & quadrature_points_coordinates,
			const ElementTypeMapArray<UInt> * filter_elements) const {


			const Array<Real> & nodes_coordinates = mesh.getNodes();

			interpolateOnIntegrationPoints(nodes_coordinates, quadrature_points_coordinates, filter_elements);
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::computeIntegrationPointsCoordinates(Array<Real> & quadrature_points_coordinates,
			const ElementType & type,
			const GhostType & ghost_type,
			const Array<UInt> & filter_elements) const {


			const Array<Real> & nodes_coordinates = mesh.getNodes();

			UInt spatial_dimension = mesh.getSpatialDimension();

			interpolateOnIntegrationPoints(nodes_coordinates, quadrature_points_coordinates, spatial_dimension,
			type, ghost_type, filter_elements);
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::initElementalFieldInterpolationFromIntegrationPoints(const ElementTypeMapArray<Real> & interpolation_points_coordinates,
			ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
			ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
			const ElementTypeMapArray<UInt> * element_filter) const {

			AKANTU_DEBUG_IN();

			UInt spatial_dimension = this->mesh.getSpatialDimension();

			ElementTypeMapArray<Real> quadrature_points_coordinates("quadrature_points_coordinates_for_interpolation", getID());
			mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension);
			computeIntegrationPointsCoordinates(quadrature_points_coordinates, element_filter);
			shape_functions.initElementalFieldInterpolationFromIntegrationPoints(interpolation_points_coordinates,
			interpolation_points_coordinates_matrices,
			quad_points_coordinates_inv_matrices,
			quadrature_points_coordinates,
			element_filter);

			}
			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::interpolateElementalFieldFromIntegrationPoints(const ElementTypeMapArray<Real> & field,
			const ElementTypeMapArray<Real> & interpolation_points_coordinates,
			ElementTypeMapArray<Real> & result,
			const GhostType ghost_type,
			const ElementTypeMapArray<UInt> * element_filter) const {

			ElementTypeMapArray<Real> interpolation_points_coordinates_matrices("interpolation_points_coordinates_matrices");
			ElementTypeMapArray<Real> quad_points_coordinates_inv_matrices("quad_points_coordinates_inv_matrices");

			initElementalFieldInterpolationFromIntegrationPoints(interpolation_points_coordinates,
			interpolation_points_coordinates_matrices,
			quad_points_coordinates_inv_matrices,
			element_filter);

			interpolateElementalFieldFromIntegrationPoints(field,
			interpolation_points_coordinates_matrices,
			quad_points_coordinates_inv_matrices,
			result,
			ghost_type,
			element_filter);
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::interpolateElementalFieldFromIntegrationPoints(const ElementTypeMapArray<Real> & field,
			const ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices,
			const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices,
			ElementTypeMapArray<Real> & result,
			const GhostType ghost_type,
			const ElementTypeMapArray<UInt> * element_filter) const {

			shape_functions.interpolateElementalFieldFromIntegrationPoints(field,
			interpolation_points_coordinates_matrices,
			quad_points_coordinates_inv_matrices,
			result,
			ghost_type,
			element_filter);
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct InterpolateHelper {
			template <class S>
			static void call(const S & shape_functions,
			const Vector<Real> & real_coords,
			UInt elem,
			const Matrix<Real> & nodal_values,
			Vector<Real> & interpolated,
			const ElementType & type,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define INTERPOLATE(type) \
			shape_functions.template interpolate<type>(real_coords, \
			element, \
			nodal_values, \
			interpolated, \
			ghost_type);

			#define AKANTU_SPECIALIZE_INTERPOLATE_HELPER(kind) \
			template<> \
			struct InterpolateHelper<kind> { \
			template <class S> \
			static void call(const S & shape_functions, \
			const Vector<Real> & real_coords, \
			UInt element, \
			const Matrix<Real> & nodal_values, \
			Vector<Real> & interpolated, \
			const ElementType & type, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INTERPOLATE_HELPER, \
			AKANTU_FE_ENGINE_LIST_INTERPOLATE)

			#undef AKANTU_SPECIALIZE_INTERPOLATE_HELPER
			#undef INTERPOLATE

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::interpolate(const Vector<Real> & real_coords,
			const Matrix<Real> & nodal_values,
			Vector<Real> & interpolated,
			const Element & element) const{

			AKANTU_DEBUG_IN();

			InterpolateHelper<kind>::call(shape_functions, real_coords, element.element, nodal_values, interpolated, element.type, element.ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::computeNormalsOnIntegrationPoints(const GhostType & ghost_type) {
			AKANTU_DEBUG_IN();

			computeNormalsOnIntegrationPoints(mesh.getNodes(),
			ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::computeNormalsOnIntegrationPoints(const Array<Real> & field,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_IN();

			// Real * coord = mesh.getNodes().storage();
			UInt spatial_dimension = mesh.getSpatialDimension();

			//allocate the normal arrays
			Mesh::type_iterator it = mesh.firstType(element_dimension, ghost_type, kind);
			Mesh::type_iterator end = mesh.lastType(element_dimension, ghost_type, kind);
			for(; it != end; ++it) {
			ElementType type = *it;
			UInt size = mesh.getNbElement(type, ghost_type);
			if(normals_on_integration_points.exists(type, ghost_type)) {
			normals_on_integration_points(type, ghost_type).resize(size);
			} else {
			normals_on_integration_points.alloc(size, spatial_dimension, type, ghost_type);
			}
			}

			//loop over the type to build the normals
			it = mesh.firstType(element_dimension, ghost_type, kind);
			for(; it != end; ++it) {
			Array<Real> & normals_on_quad = normals_on_integration_points(*it, ghost_type);
			computeNormalsOnIntegrationPoints(field, normals_on_quad, *it, ghost_type);
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct ComputeNormalsOnIntegrationPoints {
			template <template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind k>
			static void call(const FEEngineTemplate<I, S, k> & fem,
			const Array<Real> & field,
			Array<Real> & normal,
			const ElementType & type,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define COMPUTE_NORMALS_ON_INTEGRATION_POINTS(type) \
			fem.template computeNormalsOnIntegrationPoints<type>(field, normal, ghost_type);

			#define AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS(kind) \
			template<> \
			struct ComputeNormalsOnIntegrationPoints<kind> { \
			template <template <ElementKind> class I, \
			template <ElementKind> class S, \
			ElementKind k> \
			static void call(const FEEngineTemplate<I, S, k> & fem, \
			const Array<Real> & field, \
			Array<Real> & normal, \
			const ElementType & type, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_INTEGRATION_POINTS, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS, \
			AKANTU_FE_ENGINE_LIST_COMPUTE_NORMALS_ON_INTEGRATION_POINTS)

			#undef AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS
			#undef COMPUTE_NORMALS_ON_INTEGRATION_POINTS

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::computeNormalsOnIntegrationPoints(const Array<Real> & field,
			Array<Real> & normal,
			const ElementType & type,
			const GhostType & ghost_type) const {
			ComputeNormalsOnIntegrationPoints<kind>::call(*this,
			field,
			normal,
			type,
			ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			template<ElementType type>
			void FEEngineTemplate<I, S, kind>::computeNormalsOnIntegrationPoints(const Array<Real> & field,
			Array<Real> & normal,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = mesh.getSpatialDimension();
			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			UInt nb_points = getNbIntegrationPoints(type, ghost_type);

			UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize();
			normal.resize(nb_element * nb_points);
			Array<Real>::matrix_iterator normals_on_quad = normal.begin_reinterpret(spatial_dimension,
			nb_points,
			nb_element);
			Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element);
			FEEngine::extractNodalToElementField(mesh, field, f_el, type, ghost_type);

			const Matrix<Real> & quads =
			integrator. template getQuadraturePoints<type>(ghost_type);

			Array<Real>::matrix_iterator f_it = f_el.begin(spatial_dimension, nb_nodes_per_element);

			for (UInt elem = 0; elem < nb_element; ++elem) {
			ElementClass<type>::computeNormalsOnNaturalCoordinates(quads,
			*f_it,
			*normals_on_quad);
			++normals_on_quad;
			++f_it;
			}

			AKANTU_DEBUG_OUT();
			}


			/* -------------------------------------------------------------------------- */
			/* Matrix lumping functions */
			/* -------------------------------------------------------------------------- */

			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct AssembleLumpedTemplateHelper { };

			#define ASSEMBLE_LUMPED(type) \
			fem.template assembleLumpedTemplate<type>(field_1, nb_degree_of_freedom,lumped, equation_number,ghost_type)

			#define AKANTU_SPECIALIZE_ASSEMBLE_HELPER(kind) \
			template<> \
			struct AssembleLumpedTemplateHelper<kind> { \
			template <template <ElementKind> class I, \
			template <ElementKind> class S, \
			ElementKind k> \
			static void call(const FEEngineTemplate<I, S, k> & fem, \
			const Array<Real> & field_1, \
			UInt nb_degree_of_freedom, \
			Array<Real> & lumped, \
			const Array<Int> & equation_number, \
			ElementType type, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_LUMPED, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_ASSEMBLE_HELPER)

			#undef AKANTU_SPECIALIZE_ASSEMBLE_HELPER
			#undef ASSEMBLE_LUMPED

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::assembleFieldLumped(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			ElementType type,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			AssembleLumpedTemplateHelper<kind>::call(*this, field_1,
			nb_degree_of_freedom,
			lumped,
			equation_number,
			type,
			ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct AssembleFieldMatrixHelper { };

			#define ASSEMBLE_MATRIX(type) \
			fem.template assembleFieldMatrix<type>(field_1, nb_degree_of_freedom, \
			matrix, \
			ghost_type)

			#define AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER(kind) \
			template<> \
			struct AssembleFieldMatrixHelper<kind> { \
			template <template <ElementKind> class I, \
			template <ElementKind> class S, \
			ElementKind k> \
			static void call(const FEEngineTemplate<I, S, k> & fem, \
			const Array<Real> & field_1, \
			UInt nb_degree_of_freedom, \
			SparseMatrix & matrix, \
			ElementType type, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(ASSEMBLE_MATRIX, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER)

			#undef AKANTU_SPECIALIZE_ASSEMBLE_FIELD_MATRIX_HELPER
			#undef ASSEMBLE_MATRIX

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::assembleFieldMatrix(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			SparseMatrix & matrix,
			ElementType type,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();
			AssembleFieldMatrixHelper<kind>::call(*this, field_1, nb_degree_of_freedom, matrix, type, ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			template <ElementType type>
			void FEEngineTemplate<I, S, kind>::assembleLumpedTemplate(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			this->template assembleLumpedRowSum<type>(field_1, nb_degree_of_freedom,lumped, equation_number,ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			/**
			* @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$
			*/
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			template <ElementType type>
			void FEEngineTemplate<I, S, kind>::assembleLumpedRowSum(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			UInt shapes_size = ElementClass<type>::getShapeSize();

			Array<Real> * field_times_shapes = new Array<Real>(0, shapes_size * nb_degree_of_freedom);
			Array<Real> * field = new Array<Real>(field_1.getSize(), nb_degree_of_freedom);

			Array<Real>::const_scalar_iterator f1_it = field_1.begin();
			Array<Real>::const_scalar_iterator f1_end = field_1.end();
			Array<Real>::vector_iterator f_it = field->begin(nb_degree_of_freedom);

			for(;f1_it != f1_end; ++f1_it, ++f_it) {
			f_it->set(*f1_it);
			}

			shape_functions.template fieldTimesShapes<type>(field, field_times_shapes, ghost_type);
			delete field;

			UInt nb_element = mesh.getNbElement(type, ghost_type);
			Array<Real> * int_field_times_shapes = new Array<Real>(nb_element, shapes_size * nb_degree_of_freedom,
			"inte_rho_x_shapes");
			integrator.template integrate<type>(field_times_shapes, int_field_times_shapes,
			nb_degree_of_freedom * shapes_size, ghost_type, empty_filter);
			delete field_times_shapes;

			assembleArray(*int_field_times_shapes, lumped, equation_number,nb_degree_of_freedom, type, ghost_type);
			delete int_field_times_shapes;

			AKANTU_DEBUG_OUT();
			}


			/* -------------------------------------------------------------------------- */
			/**
			* @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
			*/
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			template <ElementType type>
			void FEEngineTemplate<I, S, kind>::assembleLumpedDiagonalScaling(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			const ElementType & type_p1 = ElementClass<type>::getP1ElementType();
			UInt nb_nodes_per_element_p1 = Mesh::getNbNodesPerElement(type_p1);
			UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			UInt nb_quadrature_points = integrator.template getQuadraturePoints<type>(ghost_type).cols();

			UInt nb_element = field_1.getSize() / nb_quadrature_points;
			Vector<Real> nodal_factor(nb_nodes_per_element);

			#define ASSIGN_WEIGHT_TO_NODES(corner,mid) \
			{ \
			for (UInt n = 0; n < nb_nodes_per_element_p1; n++) \
			nodal_factor(n) = corner; \
			for (UInt n = nb_nodes_per_element_p1; \
			n < nb_nodes_per_element; n++) \
			nodal_factor(n) = mid; \
			}


			if (type == _triangle_6 ) ASSIGN_WEIGHT_TO_NODES(1./12., 1./4.);
			if (type == _tetrahedron_10) ASSIGN_WEIGHT_TO_NODES(1./32., 7./48.);
			if (type == _quadrangle_8 ) ASSIGN_WEIGHT_TO_NODES(3./76., 16./76.); /** coeff. derived by scaling the diagonal terms of the corresponding
			* consistent mass computed with 3x3 gauss points;
			* coeff. are (1./36., 8./36.) for 2x2 gauss points */
			if (type == _hexahedron_20 ) ASSIGN_WEIGHT_TO_NODES(7./248., 16./248.); /** coeff. derived by scaling the diagonal terms of the corresponding
			* consistent mass computed with 3x3x3 gauss points;
			* coeff. are (1./40., 1./15.) for 2x2x2 gauss points */
			if (type == _pentahedron_15) {
			// coefficients derived by scaling the diagonal terms of the corresponding consistent mass computed with 8 gauss points;
			for (UInt n = 0; n < nb_nodes_per_element_p1; n++)
			nodal_factor(n) = 51./2358.;

			Real mid_triangle = 192./2358.;
			Real mid_quadrangle = 300./2358.;

			nodal_factor(6) = mid_triangle;
			nodal_factor(7) = mid_triangle;
			nodal_factor(8) = mid_triangle;
			nodal_factor(9) = mid_quadrangle;
			nodal_factor(10) = mid_quadrangle;
			nodal_factor(11) = mid_quadrangle;
			nodal_factor(12) = mid_triangle;
			nodal_factor(13) = mid_triangle;
			nodal_factor(14) = mid_triangle;
			}

			if (nb_element == 0) {
			AKANTU_DEBUG_OUT();
			return;
			}

			#undef ASSIGN_WEIGHT_TO_NODES

			/// compute @f$ \int \rho dV = \rho V @f$ for each element
			Array<Real> * int_field_1 = new Array<Real>(field_1.getSize(), 1,
			"inte_rho_x_1");
			integrator.template integrate<type>(field_1, *int_field_1, 1, ghost_type, empty_filter);

			/// distribute the mass of the element to the nodes
			Array<Real> * lumped_per_node = new Array<Real>(nb_element, nb_degree_of_freedom * nb_nodes_per_element, "mass_per_node");
			Array<Real>::const_scalar_iterator int_field_1_it = int_field_1->begin();
			Array<Real>::matrix_iterator lumped_per_node_it
			= lumped_per_node->begin(nb_degree_of_freedom, nb_nodes_per_element);

			for (UInt e = 0; e < nb_element; ++e) {
			for (UInt n = 0; n < nb_nodes_per_element; ++n) {
			Real lmass = int_field_1_it nodal_factor(n);
			Vector<Real> l = (*lumped_per_node_it)(n);
			l.set(lmass);
			}
			++int_field_1_it;
			++lumped_per_node_it;
			}
			delete int_field_1;

			// lumped_per_node->extendComponentsInterlaced(nb_degree_of_freedom,1);
			assembleArray(*lumped_per_node, lumped, equation_number, nb_degree_of_freedom, type, ghost_type);
			delete lumped_per_node;

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$
			*/
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			template <ElementType type>
			void FEEngineTemplate<I, S, kind>::assembleFieldMatrix(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			SparseMatrix & matrix,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			UInt vect_size = field_1.getSize();
			UInt shapes_size = ElementClass<type>::getShapeSize();
			UInt lmat_size = nb_degree_of_freedom * shapes_size;

			const Array<Real> & shapes = shape_functions.getShapes(type,ghost_type);
			Array<Real> * modified_shapes = new Array<Real>(vect_size, lmat_size * nb_degree_of_freedom);
			modified_shapes->clear();
			Array<Real> * local_mat = new Array<Real>(vect_size, lmat_size * lmat_size);

			Array<Real>::matrix_iterator mshapes_it = modified_shapes->begin(lmat_size, nb_degree_of_freedom);
			Array<Real>::const_vector_iterator shapes_it = shapes.begin(shapes_size);

			for(UInt q = 0; q < vect_size; ++q, ++mshapes_it, ++shapes_it) {
			for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
			for (UInt s = 0; s < shapes_size; ++s) {
			(mshapes_it)(snb_degree_of_freedom + d, d) = (*shapes_it)(s);
			}
			}
			}

			mshapes_it = modified_shapes->begin(lmat_size, nb_degree_of_freedom);
			Array<Real>::matrix_iterator lmat = local_mat->begin(lmat_size, lmat_size);
			Array<Real>::const_scalar_iterator field_val = field_1.begin();

			for(UInt q = 0; q < vect_size; ++q, ++lmat, ++mshapes_it, ++field_val) {
			(lmat).mul<false, true>(mshapes_it, mshapes_it, field_val);
			}

			delete modified_shapes;

			UInt nb_element = mesh.getNbElement(type, ghost_type);
			Array<Real> * int_field_times_shapes = new Array<Real>(nb_element, lmat_size * lmat_size,
			"inte_rho_x_shapes");
			integrator.template integrate<type>(local_mat, int_field_times_shapes,
			lmat_size * lmat_size, ghost_type, empty_filter);
			delete local_mat;

			assembleMatrix(*int_field_times_shapes, matrix, nb_degree_of_freedom, type, ghost_type);
			delete int_field_times_shapes;

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct InverseMapHelper {
			template <class S>
			static void call(const S & shape_functions,
			const Vector<Real> & real_coords,
			UInt element,
			const ElementType & type,
			Vector<Real> & natural_coords,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define INVERSE_MAP(type) \
			shape_functions.template inverseMap<type>(real_coords, element, natural_coords, ghost_type); \

			#define AKANTU_SPECIALIZE_INVERSE_MAP_HELPER(kind) \
			template<> \
			struct InverseMapHelper<kind> { \
			template <class S> \
			static void call(const S & shape_functions, \
			const Vector<Real> & real_coords, \
			UInt element, \
			const ElementType & type, \
			Vector<Real> & natural_coords, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(INVERSE_MAP, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INVERSE_MAP_HELPER, \
			AKANTU_FE_ENGINE_LIST_INVERSE_MAP)

			#undef AKANTU_SPECIALIZE_INVERSE_MAP_HELPER
			#undef INVERSE_MAP

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::inverseMap(const Vector<Real> & real_coords,
			UInt element,
			const ElementType & type,
			Vector<Real> & natural_coords,
			const GhostType & ghost_type) const{

			AKANTU_DEBUG_IN();

			InverseMapHelper<kind>::call(shape_functions, real_coords, element, type, natural_coords, ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct ContainsHelper {
			template <class S>
			static void call(const S & shape_functions,
			const Vector<Real> & real_coords,
			UInt element,
			const ElementType & type,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define CONTAINS(type) \
			contain = shape_functions.template contains<type>(real_coords, element, ghost_type); \

			#define AKANTU_SPECIALIZE_CONTAINS_HELPER(kind) \
			template<> \
			struct ContainsHelper<kind> { \
			template <template <ElementKind> class S, \
			ElementKind k> \
			static bool call(const S<k> & shape_functions, \
			const Vector<Real> & real_coords, \
			UInt element, \
			const ElementType & type, \
			const GhostType & ghost_type) { \
			bool contain = false; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(CONTAINS, kind); \
			return contain; \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_CONTAINS_HELPER, \
			AKANTU_FE_ENGINE_LIST_CONTAINS)

			#undef AKANTU_SPECIALIZE_CONTAINS_HELPER
			#undef CONTAINS

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline bool FEEngineTemplate<I, S, kind>::contains(const Vector<Real> & real_coords,
			UInt element,
			const ElementType & type,
			const GhostType & ghost_type) const{
			return ContainsHelper<kind>::call(shape_functions, real_coords, element, type, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct ComputeShapesHelper {
			template <class S>
			static void call(const S & shape_functions,
			const Vector<Real> & real_coords,
			UInt element,
			const ElementType type,
			Vector<Real> & shapes,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define COMPUTE_SHAPES(type) \
			shape_functions.template computeShapes<type>(real_coords, \
			element, \
			shapes, \
			ghost_type);

			#define AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER(kind) \
			template<> \
			struct ComputeShapesHelper<kind> { \
			template <class S> \
			static void call(const S & shape_functions, \
			const Vector<Real> & real_coords, \
			UInt element, \
			const ElementType type, \
			Vector<Real> & shapes, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPES, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER, \
			AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES)

			#undef AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER
			#undef COMPUTE_SHAPES

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::computeShapes(const Vector<Real> & real_coords,
			UInt element,
			const ElementType & type,
			Vector<Real> & shapes,
			const GhostType & ghost_type) const{

			AKANTU_DEBUG_IN();

			ComputeShapesHelper<kind>::call(shape_functions, real_coords, element, type, shapes, ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct ComputeShapeDerivativesHelper {
			template <class S>
			static void call(const S & shape_functions,
			const Vector<Real> & real_coords,
			UInt element,
			const ElementType type,
			Matrix<Real> & shape_derivatives,
			const GhostType & ghost_type) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define COMPUTE_SHAPE_DERIVATIVES(type) \
			Matrix<Real> coords_mat(real_coords.storage(), shape_derivatives.rows(), 1); \
			Tensor3<Real> shapesd_tensor(shape_derivatives.storage(), \
			shape_derivatives.rows(), \
			shape_derivatives.cols(), 1); \
			shape_functions.template computeShapeDerivatives<type>(coords_mat, \
			element, \
			shapesd_tensor, \
			ghost_type);

			#define AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER(kind) \
			template<> \
			struct ComputeShapeDerivativesHelper<kind> { \
			template <class S> \
			static void call(const S & shape_functions, \
			const Vector<Real> & real_coords, \
			UInt element, \
			const ElementType type, \
			Matrix<Real> & shape_derivatives, \
			const GhostType & ghost_type) { \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPE_DERIVATIVES, kind); \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER, \
			AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES_DERIVATIVES)

			#undef AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER
			#undef COMPUTE_SHAPE_DERIVATIVES

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline void FEEngineTemplate<I, S, kind>::computeShapeDerivatives(const Vector<Real> & real_coords,
			UInt element,
			const ElementType & type,
			Matrix<Real> & shape_derivatives,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			ComputeShapeDerivativesHelper<kind>::call(shape_functions, real_coords, element, type, shape_derivatives, ghost_type);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct GetNbIntegrationPointsHelper { };

			#define GET_NB_INTEGRATION_POINTS(type) \
			nb_quad_points = \
			integrator. template getQuadraturePoints<type>(ghost_type).cols();

			#define AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER(kind) \
			template<> \
			struct GetNbIntegrationPointsHelper<kind> { \
			template <template <ElementKind> class I, \
			ElementKind k> \
			static UInt call(const I<k> & integrator, \
			const ElementType type, \
			const GhostType & ghost_type) { \
			UInt nb_quad_points = 0; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_NB_INTEGRATION_POINTS, kind); \
			return nb_quad_points; \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER)

			#undef AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER
			#undef GET_NB_INTEGRATION

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline UInt FEEngineTemplate<I, S, kind>::getNbIntegrationPoints(const ElementType & type,
			const GhostType & ghost_type) const {
			return GetNbIntegrationPointsHelper<kind>::call(integrator, type, ghost_type);
			}


			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct GetShapesHelper { };

			#define GET_SHAPES(type) \
			ret = &(shape_functions.getShapes(type, ghost_type));

			#define AKANTU_SPECIALIZE_GET_SHAPES_HELPER(kind) \
			template<> \
			struct GetShapesHelper<kind> { \
			template <class S> \
			static const Array<Real> & call(const S& shape_functions, \
			const ElementType type, \
			const GhostType & ghost_type) { \
			const Array<Real> * ret = NULL; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES, kind); \
			return *ret; \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_SHAPES_HELPER)

			#undef AKANTU_SPECIALIZE_GET_SHAPES_HELPER
			#undef GET_SHAPES

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline const Array<Real> & FEEngineTemplate<I, S, kind>::getShapes(const ElementType & type,
			const GhostType & ghost_type,
			__attribute__((unused)) UInt id) const {
			return GetShapesHelper<kind>::call(shape_functions, type, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct GetShapesDerivativesHelper {
			template <template <ElementKind> class S,
			ElementKind k>
			static const Array<Real> & call(const S<k> & shape_functions,
			const ElementType & type,
			const GhostType & ghost_type,
			UInt id) {
			AKANTU_DEBUG_TO_IMPLEMENT();
			}
			};

			#define GET_SHAPES_DERIVATIVES(type) \
			ret = &(shape_functions.getShapesDerivatives(type, ghost_type));

			#define AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER(kind) \
			template<> \
			struct GetShapesDerivativesHelper<kind> { \
			template <template <ElementKind> class S, \
			ElementKind k> \
			static const Array<Real> & call(const S<k> & shape_functions, \
			const ElementType type, \
			const GhostType & ghost_type, \
			UInt id) { \
			const Array<Real> * ret = NULL; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES_DERIVATIVES, kind); \
			return *ret; \
			} \
			};

			AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER, \
			AKANTU_FE_ENGINE_LIST_GET_SHAPES_DERIVATIVES)

			#undef AKANTU_SPECIALIZE_GET_SHAPE_DERIVATIVES_HELPER
			#undef GET_SHAPES_DERIVATIVES

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline const Array<Real> & FEEngineTemplate<I, S, kind>::getShapesDerivatives(const ElementType & type,
			const GhostType & ghost_type,
			__attribute__((unused)) UInt id) const {
			return GetShapesDerivativesHelper<kind>::call(shape_functions, type, ghost_type, id);
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Helper class to be able to write a partial specialization on the element kind
			*/
			template<ElementKind kind>
			struct GetIntegrationPointsHelper { };

			#define GET_INTEGRATION_POINTS(type) \
			ret = &(integrator. template getQuadraturePoints<type>(ghost_type)); \

			#define AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER(kind) \
			template<> \
			struct GetIntegrationPointsHelper<kind> { \
			template <template <ElementKind> class I, \
			ElementKind k> \
			static const Matrix<Real> & call(const I<k> & integrator, \
			const ElementType type, \
			const GhostType & ghost_type) { \
			const Matrix<Real> * ret = NULL; \
			AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_INTEGRATION_POINTS, kind); \
			return *ret; \
			} \
			};

			AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER)

			#undef AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER
			#undef GET_INTEGRATION_POINTS

			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			inline const Matrix<Real> &
			FEEngineTemplate<I, S, kind>::getIntegrationPoints(const ElementType & type,
			const GhostType & ghost_type) const {
			return GetIntegrationPointsHelper<kind>::call(integrator, type, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template<template <ElementKind> class I,
			template <ElementKind> class S,
			ElementKind kind>
			void FEEngineTemplate<I, S, kind>::printself(std::ostream & stream, int indent) const {
			std::string space;
			for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);

			stream << space << "FEEngineTemplate [" << std::endl;
			stream << space << " + parent [" << std::endl;
			FEEngine::printself(stream, indent + 3);
			stream << space << " ]" << std::endl;
			stream << space << " + shape functions [" << std::endl;
			shape_functions.printself(stream, indent + 3);
			stream << space << " ]" << std::endl;
			stream << space << " + integrator [" << std::endl;
			integrator.printself(stream, indent + 3);
			stream << space << " ]" << std::endl;
			stream << space << "]" << std::endl;
			}

			/* -------------------------------------------------------------------------- */

			__END_AKANTU__

			#include "shape_lagrange.hh"
			#include "integrator_gauss.hh"

			__BEGIN_AKANTU__

			/* -------------------------------------------------------------------------- */
			template <>
			template <>
			inline void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
			assembleLumpedTemplate<_triangle_6>(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			assembleLumpedDiagonalScaling<_triangle_6>(field_1, nb_degree_of_freedom,
			lumped, equation_number, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template <>
			template <>
			inline void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
			assembleLumpedTemplate<_tetrahedron_10>(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			assembleLumpedDiagonalScaling<_tetrahedron_10>(field_1, nb_degree_of_freedom,
			lumped, equation_number, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template <>
			template <>
			inline void
			FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
			assembleLumpedTemplate<_quadrangle_8>(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			assembleLumpedDiagonalScaling<_quadrangle_8>(field_1, nb_degree_of_freedom,
			lumped, equation_number, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template <>
			template <>
			inline void
			FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
			assembleLumpedTemplate<_hexahedron_20>(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			assembleLumpedDiagonalScaling<_hexahedron_20>(field_1, nb_degree_of_freedom,
			lumped, equation_number, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template <>
			template <>
			inline void
			FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
			assembleLumpedTemplate<_pentahedron_15>(const Array<Real> & field_1,
			UInt nb_degree_of_freedom,
			Array<Real> & lumped,
			const Array<Int> & equation_number,
			const GhostType & ghost_type) const {
			assembleLumpedDiagonalScaling<_pentahedron_15>(field_1, nb_degree_of_freedom,
			lumped, equation_number, ghost_type);
			}

			/* -------------------------------------------------------------------------- */
			template<>
			template<>
			inline void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
			computeNormalsOnIntegrationPoints<_point_1>(__attribute__((unused))const Array<Real> & field,
			Array<Real> & normal,
			const GhostType & ghost_type) const {
			AKANTU_DEBUG_IN();

			AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == 1, "Mesh dimension must be 1 to compute normals on points!");
			const ElementType type = _point_1;
			UInt spatial_dimension = mesh.getSpatialDimension();
			//UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			UInt nb_points = getNbIntegrationPoints(type, ghost_type);

			UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize();
			normal.resize(nb_element * nb_points);
			Array<Real>::matrix_iterator normals_on_quad = normal.begin_reinterpret(spatial_dimension,
			nb_points,
			nb_element);
			Array< std::vector<Element> > segments = mesh.getElementToSubelement(type, ghost_type);
			Array<Real> coords = mesh.getNodes();

			const Mesh * mesh_segment;
			if (mesh.isMeshFacets())
			mesh_segment = &(mesh.getMeshParent());
			else
			mesh_segment = &mesh;

			for (UInt elem = 0; elem < nb_element; ++elem) {
			UInt nb_segment = segments(elem).size();

			AKANTU_DEBUG_ASSERT(nb_segment > 0, "Impossible to compute a normal on a point connected to 0 segments");

			Real normal_value = 1;
			if (nb_segment == 1) {
			Element segment = segments(elem)[0];
			const Array<UInt> & segment_connectivity = mesh_segment->getConnectivity(segment.type, segment.ghost_type);
			//const Vector<UInt> & segment_points = segment_connectivity.begin(Mesh::getNbNodesPerElement(segment.type))[segment.element];
			Real difference;
			if(segment_connectivity(0) == elem) {
			difference = coords(elem)-coords(segment_connectivity(1));
			} else {
			difference = coords(elem)-coords(segment_connectivity(0));
			}
			normal_value = difference/std::abs(difference);
			}

			for(UInt n(0); n < nb_points; ++n) {
			(*normals_on_quad)(0, n) = normal_value;
			}
			++normals_on_quad;
			}

			AKANTU_DEBUG_OUT();
			}

			__END_AKANTU__

rAKA/src/fe_engine/fe_engine_template_tmpl.hh54e599859086bugfixes/bug_fix_in_syn…

fe_engine_template_tmpl.hhView Options

rAKA/src/fe_engine/fe_engine_template_tmpl.hh
54e599859086bugfixes/bug_fix_in_syn…

fe_engine_template_tmpl.hh
View Options