diff --git a/src/fe_engine/shape_cohesive.hh b/src/fe_engine/shape_cohesive.hh index d214f7035..5e61e7e2f 100644 --- a/src/fe_engine/shape_cohesive.hh +++ b/src/fe_engine/shape_cohesive.hh @@ -1,169 +1,183 @@ /** * @file shape_cohesive.hh * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Oct 22 2015 * * @brief shape functions for cohesive elements description * * @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_array.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_COHESIVE_HH__ #define __AKANTU_SHAPE_COHESIVE_HH__ namespace akantu { struct CohesiveReduceFunctionMean { inline Real operator()(Real u_plus, Real u_minus) { return .5 * (u_plus + u_minus); } }; struct CohesiveReduceFunctionOpening { inline Real operator()(Real u_plus, Real u_minus) { return (u_plus - u_minus); } }; template <> class ShapeLagrange<_ek_cohesive> : public ShapeLagrangeBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ShapeLagrange(const Mesh & mesh, const ID & id = "shape_cohesive", const MemoryID & memory_id = 0); ~ShapeLagrange() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initShapeFunctions(const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type); /// extract the nodal values and store them per element template void extractNodalToElementField( const Array & nodal_f, Array & elemental_f, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// computes the shape functions derivatives for given interpolation points template void computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const GhostType & ghost_type, const Array & filter_elements = empty_filter) const; void computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const ElementType & type, const GhostType & ghost_type, const Array & filter_elements) const override; /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// pre compute all shape derivatives on the element integration points from /// natural coordinates template void precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// interpolate nodal values on the integration points template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, const GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, const GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const { interpolateOnIntegrationPoints( u, uq, nb_degree_of_freedom, ghost_type, filter_elements); } /// compute the gradient of u on the integration points in the natural /// coordinates template void gradientOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const { variationOnIntegrationPoints( u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements); } + template + void computeBtD(const Array & /*Ds*/, Array & /*BtDs*/, + GhostType /*ghost_type*/, + const Array & /*filter_elements*/) const { + AKANTU_DEBUG_TO_IMPLEMENT(); + } + + template + void computeBtDB(const Array & /*Ds*/, Array & /*BtDBs*/, + UInt /*order_d*/, GhostType /*ghost_type*/, + const Array & /*filter_elements*/) const { + AKANTU_DEBUG_TO_IMPLEMENT(); + } + /// compute the gradient of u on the integration points template void variationOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// compute the normals to the field u on integration points template void computeNormalsOnIntegrationPoints( const Array & u, Array & normals_u, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// multiply a field by shape functions template void fieldTimesShapes(__attribute__((unused)) const Array & field, __attribute__((unused)) Array & fiedl_times_shapes, __attribute__((unused)) GhostType ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } }; /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const ShapeCohesive & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "shape_cohesive_inline_impl.cc" #endif /* __AKANTU_SHAPE_COHESIVE_HH__ */ diff --git a/src/fe_engine/shape_lagrange_inline_impl.cc b/src/fe_engine/shape_lagrange_inline_impl.cc index a4b2cc5ea..f2cf401c2 100644 --- a/src/fe_engine/shape_lagrange_inline_impl.cc +++ b/src/fe_engine/shape_lagrange_inline_impl.cc @@ -1,528 +1,526 @@ /** * @file shape_lagrange_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Wed Oct 27 2010 * @date last modification: Thu Oct 15 2015 * * @brief ShapeLagrange inline implementation * * @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_iterators.hh" #include "aka_voigthelper.hh" #include "fe_engine.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ #define __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ if (ElementClass::getNaturalSpaceDimension() == \ mesh.getSpatialDimension() || \ kind != _ek_regular) \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); template inline void ShapeLagrange::initShapeFunctions( const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type) { AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } #undef INIT_SHAPE_FUNCTIONS /* -------------------------------------------------------------------------- */ template template inline void ShapeLagrange::computeShapeDerivativesOnCPointsByElement( const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const { AKANTU_DEBUG_IN(); // compute dnds Tensor3 dnds(node_coords.rows(), node_coords.cols(), natural_coords.cols()); ElementClass::computeDNDS(natural_coords, dnds); // compute jacobian Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); // compute dndx ElementClass::computeShapeDerivatives(J, dnds, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); ElementClass::inverseMap(real_coords, nodes_coord, natural_coords); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template bool ShapeLagrange::contains(const Vector & real_coords, UInt elem, const GhostType & ghost_type) const { UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); return ElementClass::contains(natural_coords); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, Vector & interpolated, const GhostType & ghost_type) const { UInt nb_shapes = ElementClass::getShapeSize(); Vector shapes(nb_shapes); computeShapes(real_coords, elem, shapes, ghost_type); ElementClass::interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); ElementClass::computeShapes(natural_coords, shapes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapeDerivatives( const Matrix & real_coords, UInt elem, Tensor3 & shapesd, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_points = real_coords.cols(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == shapesd.size(0) && nb_nodes_per_element == shapesd.size(1), "Shape size doesn't match"); AKANTU_DEBUG_ASSERT(nb_points == shapesd.size(2), "Number of points doesn't match shapes size"); Matrix natural_coords(spatial_dimension, nb_points); // Creates the matrix of natural coordinates for (UInt i = 0; i < nb_points; i++) { Vector real_point = real_coords(i); Vector natural_point = natural_coords(i); inverseMap(real_point, elem, natural_point, ghost_type); } UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); computeShapeDerivativesOnCPointsByElement(nodes_coord, natural_coords, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template ShapeLagrange::ShapeLagrange(const Mesh & mesh, const ID & id, const MemoryID & memory_id) : ShapeLagrangeBase(mesh, kind, id, memory_id) {} /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd, "The shapes_derivatives array does not have the correct " << "number of component"); shape_derivatives.resize(nb_element * nb_points); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type, filter_elements); Real * shapesd_val = shape_derivatives.storage(); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); if (filter_elements != empty_filter) nb_element = filter_elements.size(); for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) { if (filter_elements != empty_filter) shapesd_val = shape_derivatives.storage() + filter_elements(elem) * size_of_shapesd * nb_points; Matrix & X = *x_it; Tensor3 B(shapesd_val, spatial_dimension, nb_nodes_per_element, nb_points); computeShapeDerivativesOnCPointsByElement(X, integration_points, B); if (filter_elements == empty_filter) shapesd_val += size_of_shapesd * nb_points; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const ElementType & type, const GhostType & ghost_type, const Array & filter_elements) const { #define AKANTU_COMPUTE_SHAPES(type) \ computeShapeDerivativesOnIntegrationPoints( \ nodes, integration_points, shape_derivatives, ghost_type, \ filter_elements); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES); #undef AKANTU_COMPUTE_SHAPES } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::precomputeShapesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapes = ElementClass::getShapeSize(); Array & shapes_tmp = shapes.alloc(0, size_of_shapes, itp_type, ghost_type); this->computeShapesOnIntegrationPoints(nodes, natural_coords, shapes_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); Array & shapes_derivatives_tmp = shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); this->computeShapeDerivativesOnIntegrationPoints( nodes, natural_coords, shapes_derivatives_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, const Array & shapes, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->interpolateElementalFieldOnIntegrationPoints( u_el, out_uq, ghost_type, shapes, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type), "No shapes for the type " << shapes.printType(itp_type, ghost_type)); this->interpolateOnIntegrationPoints(in_u, out_uq, nb_degree_of_freedom, shapes(itp_type, ghost_type), ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT( shapes_derivatives.exists(itp_type, ghost_type), "No shapes derivatives for the type " << shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->gradientElementalFieldOnIntegrationPoints( u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeBtD( const Array & Ds, Array & BtDs, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); auto spatial_dimension = mesh.getSpatialDimension(); - auto size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, spatial_dimension, nb_nodes_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, spatial_dimension, nb_nodes_per_element); B_it = view.begin(); B_end = view.end(); } for (auto && values : zip(range(B_it, B_end), make_view(Ds, Ds.getNbComponent() / spatial_dimension, spatial_dimension), make_view(BtDs, BtDs.getNbComponent() / nb_nodes_per_element, nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D = std::get<2>(values); // transposed due to the storage layout of B Bt_D.template mul(D, B); } } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::computeBtDB( const Array & Ds, Array & BtDBs, UInt order_d, GhostType ghost_type, const Array & filter_elements) const { auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); constexpr auto dim = ElementClass::getSpatialDimension(); - auto size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, dim, nb_nodes_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, dim, nb_nodes_per_element); B_it = view.begin(); B_end = view.end(); } if (order_d == 4) { UInt tangent_size = VoigtHelper::size; Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); for (auto && values : zip(range(B_it, B_end), make_view(Ds, tangent_size, tangent_size), make_view(BtDBs, dim * nb_nodes_per_element, dim * nb_nodes_per_element))) { const auto & Bfull = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D_B = std::get<2>(values); VoigtHelper::transferBMatrixToSymVoigtBMatrix(Bfull, B, nb_nodes_per_element); Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } } else if (order_d == 2) { Matrix Bt_D(nb_nodes_per_element, dim); for (auto && values : zip(range(B_it, B_end), make_view(Ds, dim, dim), make_view(BtDBs, nb_nodes_per_element, nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); auto & Bt_D_B = std::get<2>(values); Bt_D.template mul(B, D); Bt_D_B.template mul(Bt_D, B); } } } template <> template <> inline void ShapeLagrange<_ek_regular>::computeBtDB<_point_1>( const Array & /*Ds*/, Array & /*BtDBs*/, UInt /*order_d*/, GhostType /*ghost_type*/, const Array & /*filter_elements*/) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template template void ShapeLagrange::fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const { AKANTU_DEBUG_IN(); field_times_shapes.resize(field.size()); UInt size_of_shapes = ElementClass::getShapeSize(); InterpolationType itp_type = ElementClassProperty::interpolation_type; UInt nb_degree_of_freedom = field.getNbComponent(); const auto & shape = shapes(itp_type, ghost_type); auto field_it = field.begin(nb_degree_of_freedom, 1); auto shapes_it = shape.begin(1, size_of_shapes); auto it = field_times_shapes.begin(nb_degree_of_freedom, size_of_shapes); auto end = field_times_shapes.end(nb_degree_of_freedom, size_of_shapes); for (; it != end; ++it, ++field_it, ++shapes_it) { it->template mul(*field_it, *shapes_it); } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ */