Diffusion akantu (bugfixes/27-bug-in-the-python-interface) Browse

27-bug-in-the…

shape_structural_inline_impl.hh
View Options

			/**
			* @file shape_structural_inline_impl.hh
			*
			* @author Fabian Barras <fabian.barras@epfl.ch>
			* @author Lucas Frerot <lucas.frerot@epfl.ch>
			* @author Nicolas Richart <nicolas.richart@epfl.ch>
			*
			* @date creation: Wed Oct 11 2017
			* @date last modification: Wed Feb 21 2018
			*
			* @brief ShapeStructural inline implementation
			*
			*
			* Copyright (©) 2016-2018 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 "mesh_iterators.hh"
			#include "shape_structural.hh"
			/* -------------------------------------------------------------------------- */

			#ifndef AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_
			#define AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_

			namespace akantu {

			namespace {
			/// Extract nodal coordinates per elements
			template <ElementType type>
			std::unique_ptr<Array<Real>> getNodesPerElement(const Mesh & mesh,
			const Array<Real> & nodes,
			GhostType ghost_type) {
			const auto dim = ElementClass<type>::getSpatialDimension();
			const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);

			auto nodes_per_element =
			std::make_unique<Array<Real>>(0, dim * nb_nodes_per_element);
			FEEngine::extractNodalToElementField(mesh, nodes, *nodes_per_element, type,
			ghost_type);
			return nodes_per_element;
			}
			} // namespace

			template <ElementKind kind>
			inline void ShapeStructural<kind>::initShapeFunctions(
			const Array<Real> & /* unused /, const Matrix<Real> & / unused */,
			ElementType /* unused /, GhostType / unused */) {
			AKANTU_TO_IMPLEMENT();
			}

			/* -------------------------------------------------------------------------- */
			#define INIT_SHAPE_FUNCTIONS(type) \
			setIntegrationPointsByType<type>(integration_points, ghost_type); \
			precomputeRotationMatrices<type>(nodes, ghost_type); \
			precomputeShapesOnIntegrationPoints<type>(nodes, ghost_type); \
			precomputeShapeDerivativesOnIntegrationPoints<type>(nodes, ghost_type);

			template <>
			inline void ShapeStructural<_ek_structural>::initShapeFunctions(
			const Array<Real> & nodes, const Matrix<Real> & integration_points,
			ElementType type, GhostType ghost_type) {
			AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
			}

			#undef INIT_SHAPE_FUNCTIONS

			/* -------------------------------------------------------------------------- */
			template <ElementKind kind>
			template <ElementType type>
			void ShapeStructural<kind>::computeShapesOnIntegrationPointsInternal(
			const Array<Real> & nodes, const Matrix<Real> & integration_points,
			Array<Real> & shapes, GhostType ghost_type,
			const Array<UInt> & filter_elements, bool mass) const {

			auto nb_points = integration_points.cols();
			auto nb_element = mesh.getConnectivity(type, ghost_type).size();
			auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();

			shapes.resize(nb_element * nb_points);

			auto nb_dofs = ElementClass<type>::getNbDegreeOfFreedom();
			auto nb_rows = nb_dofs;
			if (mass) {
			nb_rows = ElementClass<type>::getNbStressComponents();
			}

			#if !defined(AKANTU_NDEBUG)
			UInt size_of_shapes = nb_rows * nb_dofs * nb_nodes_per_element;
			AKANTU_DEBUG_ASSERT(shapes.getNbComponent() == size_of_shapes,
			"The shapes array does not have the correct "
			<< "number of component");
			#endif


			auto shapes_it = shapes.begin_reinterpret(
			nb_rows, ElementClass<type>::getNbNodesPerInterpolationElement() * nb_dofs,
			nb_points, nb_element);
			auto shapes_begin = shapes_it;
			if (filter_elements != empty_filter) {
			nb_element = filter_elements.size();
			}

			auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);
			auto nodes_it = nodes_per_element->begin(mesh.getSpatialDimension(),
			Mesh::getNbNodesPerElement(type));
			auto nodes_begin = nodes_it;
			auto rot_matrix_it =
			make_view(rotation_matrices(type, ghost_type), nb_dofs, nb_dofs).begin();
			auto rot_matrix_begin = rot_matrix_it;

			for (UInt elem = 0; elem < nb_element; ++elem) {
			if (filter_elements != empty_filter) {
			shapes_it = shapes_begin + filter_elements(elem);
			nodes_it = nodes_begin + filter_elements(elem);
			rot_matrix_it = rot_matrix_begin + filter_elements(elem);
			}

			Tensor3<Real> & N = *shapes_it;
			auto & real_coord = *nodes_it;

			auto & RDOFs = *rot_matrix_it;

			Matrix<Real> T(N.size(1), N.size(1), 0);

			for (UInt i = 0; i < nb_nodes_per_element; ++i) {
			T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
			}

			if (not mass) {
			ElementClass<type>::computeShapes(integration_points, real_coord, T, N);
			} else {
			ElementClass<type>::computeShapesMass(integration_points, real_coord, T,
			N);
			}
			if (filter_elements == empty_filter) {
			++shapes_it;
			++nodes_it;
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			template <ElementKind kind>
			template <ElementType type>
			void ShapeStructural<kind>::precomputeRotationMatrices(
			const Array<Real> & nodes, GhostType ghost_type) {
			AKANTU_DEBUG_IN();

			const auto spatial_dimension = mesh.getSpatialDimension();
			const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			const auto nb_element = mesh.getNbElement(type, ghost_type);
			const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();

			if (not this->rotation_matrices.exists(type, ghost_type)) {
			this->rotation_matrices.alloc(0, nb_dof * nb_dof, type, ghost_type);
			}

			auto & rot_matrices = this->rotation_matrices(type, ghost_type);
			rot_matrices.resize(nb_element);

			Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
			FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);

			bool has_extra_normal = mesh.hasData<Real>("extra_normal", type, ghost_type);
			Array<Real>::const_vector_iterator extra_normal;
			if (has_extra_normal) {
			extra_normal = mesh.getData<Real>("extra_normal", type, ghost_type)
			.begin(spatial_dimension);
			}

			for (auto && tuple :
			zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
			make_view(rot_matrices, nb_dof, nb_dof))) {
			// compute shape derivatives
			auto & X = std::get<0>(tuple);
			auto & R = std::get<1>(tuple);

			if (has_extra_normal) {
			ElementClass<type>::computeRotationMatrix(R, X, *extra_normal);
			++extra_normal;
			} else {
			ElementClass<type>::computeRotationMatrix(
			R, X, Vector<Real>(spatial_dimension));
			}
			}

			AKANTU_DEBUG_OUT();
			}

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

			template <ElementKind kind>
			template <ElementType type>
			void ShapeStructural<kind>::precomputeShapesOnIntegrationPoints(
			const Array<Real> & nodes, GhostType ghost_type) {
			AKANTU_DEBUG_IN();

			const auto & natural_coords = integration_points(type, ghost_type);
			auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			auto nb_points = integration_points(type, ghost_type).cols();
			auto nb_element = mesh.getNbElement(type, ghost_type);
			auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
			const auto dim = ElementClass<type>::getSpatialDimension();
			const auto spatial_dimension = mesh.getSpatialDimension();
			const auto natural_spatial_dimension =
			ElementClass<type>::getNaturalSpaceDimension();

			auto itp_type = FEEngine::getInterpolationType(type);
			if (not shapes.exists(itp_type, ghost_type)) {
			auto size_of_shapes = this->getShapeSize(type);
			this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type);
			}

			auto & rot_matrices = this->rotation_matrices(type, ghost_type);
			auto & shapes_ = this->shapes(itp_type, ghost_type);
			shapes_.resize(nb_element * nb_points);

			auto nodes_per_element = getNodesPerElement<type>(mesh, nodes, ghost_type);

			for (auto && tuple :
			zip(make_view(shapes_, nb_dof, nb_dof * nb_nodes_per_element, nb_points),
			make_view(*nodes_per_element, dim, nb_nodes_per_element),
			make_view(rot_matrices, nb_dof, nb_dof))) {
			auto & N = std::get<0>(tuple);
			auto & X = std::get<1>(tuple);
			auto & RDOFs = std::get<2>(tuple);

			Matrix<Real> T(N.size(1), N.size(1), 0);

			for (UInt i = 0; i < nb_nodes_per_element; ++i) {
			T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
			}

			auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
			// Rotate to local basis
			auto x =
			(R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);

			ElementClass<type>::computeShapes(natural_coords, x, T, N);
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <ElementKind kind>
			template <ElementType type>
			void ShapeStructural<kind>::precomputeShapeDerivativesOnIntegrationPoints(
			const Array<Real> & nodes, GhostType ghost_type) {
			AKANTU_DEBUG_IN();

			const auto & natural_coords = integration_points(type, ghost_type);
			const auto spatial_dimension = mesh.getSpatialDimension();
			const auto natural_spatial_dimension =
			ElementClass<type>::getNaturalSpaceDimension();
			const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
			const auto nb_points = natural_coords.cols();
			const auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
			const auto nb_element = mesh.getNbElement(type, ghost_type);
			const auto nb_stress_components = ElementClass<type>::getNbStressComponents();

			auto itp_type = FEEngine::getInterpolationType(type);
			if (not this->shapes_derivatives.exists(itp_type, ghost_type)) {
			auto size_of_shapesd = this->getShapeDerivativesSize(type);
			this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type);
			}

			auto & rot_matrices = this->rotation_matrices(type, ghost_type);

			Array<Real> x_el(0, spatial_dimension * nb_nodes_per_element);
			FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type);

			auto & shapesd = this->shapes_derivatives(itp_type, ghost_type);
			shapesd.resize(nb_element * nb_points);

			for (auto && tuple :
			zip(make_view(x_el, spatial_dimension, nb_nodes_per_element),
			make_view(shapesd, nb_stress_components,
			nb_nodes_per_element * nb_dof, nb_points),
			make_view(rot_matrices, nb_dof, nb_dof))) {
			// compute shape derivatives
			auto & X = std::get<0>(tuple);
			auto & B = std::get<1>(tuple);
			auto & RDOFs = std::get<2>(tuple);

			Tensor3<Real> dnds(natural_spatial_dimension,
			ElementClass<type>::interpolation_property::dnds_columns,
			B.size(2));
			ElementClass<type>::computeDNDS(natural_coords, X, dnds);

			Tensor3<Real> J(natural_spatial_dimension, natural_spatial_dimension,
			natural_coords.cols());

			// Computing the coordinates of the element in the natural space
			auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension);
			Matrix<Real> T(B.size(1), B.size(1), 0);

			for (UInt i = 0; i < nb_nodes_per_element; ++i) {
			T.block(RDOFs, i * RDOFs.rows(), i * RDOFs.rows());
			}

			// Rotate to local basis
			auto x =
			(R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element);

			ElementClass<type>::computeJMat(natural_coords, x, J);
			ElementClass<type>::computeShapeDerivatives(J, dnds, T, B);
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <ElementKind kind>
			template <ElementType type>
			void ShapeStructural<kind>::interpolateOnIntegrationPoints(
			const Array<Real> & in_u, Array<Real> & out_uq, UInt nb_dof,
			GhostType ghost_type, const Array<UInt> & filter_elements) const {
			AKANTU_DEBUG_IN();

			AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_dof,
			"The output array shape is not correct");

			auto itp_type = FEEngine::getInterpolationType(type);
			const auto & shapes_ = shapes(itp_type, ghost_type);

			auto nb_element = mesh.getNbElement(type, ghost_type);
			auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();
			auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();

			Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
			FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
			filter_elements);

			auto nb_quad_points = nb_quad_points_per_element * u_el.size();
			out_uq.resize(nb_quad_points);

			auto out_it = out_uq.begin_reinterpret(nb_dof, 1, nb_quad_points_per_element,
			u_el.size());
			auto shapes_it =
			shapes_.begin_reinterpret(nb_dof, nb_dof * nb_nodes_per_element,
			nb_quad_points_per_element, nb_element);
			auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
			nb_quad_points_per_element, u_el.size());

			for_each_element(nb_element, filter_elements, [&](auto && el) {
			auto & uq = *out_it;
			const auto & u = *u_it;
			auto N = Tensor3<Real>(shapes_it[el]);

			for (auto && q : arange(uq.size(2))) {
			auto uq_q = Matrix<Real>(uq(q));
			auto u_q = Matrix<Real>(u(q));
			auto N_q = Matrix<Real>(N(q));

			uq_q.mul<false, false>(N_q, u_q);
			}

			++out_it;
			++u_it;
			});
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <ElementKind kind>
			template <ElementType type>
			void ShapeStructural<kind>::gradientOnIntegrationPoints(
			const Array<Real> & in_u, Array<Real> & out_nablauq, UInt nb_dof,
			GhostType ghost_type, const Array<UInt> & filter_elements) const {
			AKANTU_DEBUG_IN();

			auto itp_type = FEEngine::getInterpolationType(type);
			const auto & shapesd = shapes_derivatives(itp_type, ghost_type);

			auto nb_element = mesh.getNbElement(type, ghost_type);
			auto element_dimension = ElementClass<type>::getSpatialDimension();
			auto nb_quad_points_per_element = integration_points(type, ghost_type).cols();
			auto nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement();

			Array<Real> u_el(0, nb_nodes_per_element * nb_dof);
			FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type,
			filter_elements);

			auto nb_quad_points = nb_quad_points_per_element * u_el.size();
			out_nablauq.resize(nb_quad_points);

			auto out_it = out_nablauq.begin_reinterpret(
			element_dimension, 1, nb_quad_points_per_element, u_el.size());
			auto shapesd_it = shapesd.begin_reinterpret(
			element_dimension, nb_dof * nb_nodes_per_element,
			nb_quad_points_per_element, nb_element);
			auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1,
			nb_quad_points_per_element, u_el.size());

			for_each_element(nb_element, filter_elements, [&](auto && el) {
			auto & nablau = *out_it;
			const auto & u = *u_it;
			auto B = Tensor3<Real>(shapesd_it[el]);

			for (auto && q : arange(nablau.size(2))) {
			auto nablau_q = Matrix<Real>(nablau(q));
			auto u_q = Matrix<Real>(u(q));
			auto B_q = Matrix<Real>(B(q));

			nablau_q.mul<false, false>(B_q, u_q);
			}

			++out_it;
			++u_it;
			});

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <>
			template <ElementType type>
			void ShapeStructural<_ek_structural>::computeBtD(
			const Array<Real> & Ds, Array<Real> & BtDs, GhostType ghost_type,
			const Array<UInt> & filter_elements) const {
			auto itp_type = ElementClassProperty<type>::interpolation_type;

			auto nb_stress = ElementClass<type>::getNbStressComponents();
			auto nb_dof_per_element = ElementClass<type>::getNbDegreeOfFreedom() *
			mesh.getNbNodesPerElement(type);

			const auto & shapes_derivatives =
			this->shapes_derivatives(itp_type, ghost_type);

			Array<Real> shapes_derivatives_filtered(0,
			shapes_derivatives.getNbComponent());
			auto && view = make_view(shapes_derivatives, nb_stress, nb_dof_per_element);
			auto B_it = view.begin();
			auto B_end = view.end();

			if (filter_elements != empty_filter) {
			FEEngine::filterElementalData(this->mesh, shapes_derivatives,
			shapes_derivatives_filtered, type, ghost_type,
			filter_elements);
			auto && view =
			make_view(shapes_derivatives_filtered, nb_stress, nb_dof_per_element);
			B_it = view.begin();
			B_end = view.end();
			}

			for (auto && values : zip(range(B_it, B_end), make_view(Ds, nb_stress),
			make_view(BtDs, BtDs.getNbComponent()))) {
			const auto & B = std::get<0>(values);
			const auto & D = std::get<1>(values);
			auto & Bt_D = std::get<2>(values);
			Bt_D.template mul<true>(B, D);
			}
			}

			/* -------------------------------------------------------------------------- */
			template <>
			template <ElementType type>
			void ShapeStructural<_ek_structural>::computeNtb(
			const Array<Real> & bs, Array<Real> & Ntbs, GhostType ghost_type,
			const Array<UInt> & filter_elements) const {
			auto itp_type = ElementClassProperty<type>::interpolation_type;

			auto nb_dof = ElementClass<type>::getNbDegreeOfFreedom();
			auto nb_nodes_per_element = mesh.getNbNodesPerElement(type);

			const auto & shapes = this->shapes(itp_type, ghost_type);

			Array<Real> shapes_filtered(0, shapes.getNbComponent());
			auto && view = make_view(shapes, nb_dof, nb_dof * nb_nodes_per_element);
			auto N_it = view.begin();
			auto N_end = view.end();

			if (filter_elements != empty_filter) {
			FEEngine::filterElementalData(this->mesh, shapes, shapes_filtered, type,
			ghost_type, filter_elements);
			auto && view =
			make_view(shapes_filtered, nb_dof, nb_dof * nb_nodes_per_element);
			N_it = view.begin();
			N_end = view.end();
			}

			for (auto && values : zip(range(N_it, N_end), make_view(bs, nb_dof),
			make_view(Ntbs, nb_dof * nb_nodes_per_element))) {
			const auto & N = std::get<0>(values);
			const auto & b = std::get<1>(values);
			auto & Nt_b = std::get<2>(values);
			Nt_b.template mul<true>(N, b);
			}
			}

			} // namespace akantu

			#endif /* AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_HH_ */

rAKA/src/fe_engine/shape_structural_inline_impl.hh5b3351523a74bugfixes/27-bug-in-the…

shape_structural_inline_impl.hhView Options

rAKA/src/fe_engine/shape_structural_inline_impl.hh
5b3351523a74bugfixes/27-bug-in-the…

shape_structural_inline_impl.hh
View Options