diff --git a/src/fe_engine/integrator_gauss_inline_impl.cc b/src/fe_engine/integrator_gauss_inline_impl.cc index d066b82fc..0ba03345d 100644 --- a/src/fe_engine/integrator_gauss_inline_impl.cc +++ b/src/fe_engine/integrator_gauss_inline_impl.cc @@ -1,645 +1,641 @@ /** * @file integrator_gauss_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Nov 19 2015 * * @brief inline function of gauss integrator * * @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 "fe_engine.hh" #include "mesh_iterators.hh" #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss::integrateOnElement( const Array & f, Real * intf, UInt nb_degree_of_freedom, const UInt elem, const GhostType & ghost_type) const { Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom, "The vector f do not have the good number of component."); Real * f_val = f.storage() + elem * f.getNbComponent(); Real * jac_val = jac_loc.storage() + elem * nb_quadrature_points; integrate(f_val, jac_val, intf, nb_degree_of_freedom, nb_quadrature_points); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss::integrate( const Vector & in_f, UInt index, const GhostType & ghost_type) const { const Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = GaussIntegrationElement::getNbQuadraturePoints(); AKANTU_DEBUG_ASSERT(in_f.size() == nb_quadrature_points, "The vector f do not have nb_quadrature_points entries."); Real * jac_val = jac_loc.storage() + index * nb_quadrature_points; Real intf; integrate(in_f.storage(), jac_val, &intf, 1, nb_quadrature_points); return intf; } /* -------------------------------------------------------------------------- */ template inline void IntegratorGauss::integrate( Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const { memset(inte, 0, nb_degree_of_freedom * sizeof(Real)); Real * cjac = jac; for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt dof = 0; dof < nb_degree_of_freedom; ++dof) { inte[dof] += *f * *cjac; ++f; } ++cjac; } } /* -------------------------------------------------------------------------- */ template template inline const Matrix & IntegratorGauss::getIntegrationPoints( const GhostType & ghost_type) const { AKANTU_DEBUG_ASSERT( quadrature_points.exists(type, ghost_type), "Quadrature points for type " << quadrature_points.printType(type, ghost_type) << " have not been initialized." << " Did you use 'computeQuadraturePoints' function ?"); return quadrature_points(type, ghost_type); } /* -------------------------------------------------------------------------- */ template template inline UInt IntegratorGauss::getNbIntegrationPoints( const GhostType & ghost_type) const { AKANTU_DEBUG_ASSERT( quadrature_points.exists(type, ghost_type), "Quadrature points for type " << quadrature_points.printType(type, ghost_type) << " have not been initialized." << " Did you use 'computeQuadraturePoints' function ?"); return quadrature_points(type, ghost_type).cols(); } /* -------------------------------------------------------------------------- */ template template inline Matrix IntegratorGauss::getIntegrationPoints() const { return GaussIntegrationElement::getQuadraturePoints(); } /* -------------------------------------------------------------------------- */ template template inline Vector IntegratorGauss::getIntegrationWeights() const { return GaussIntegrationElement::getWeights(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss::computeQuadraturePoints( const GhostType & ghost_type) { Matrix & quads = quadrature_points(type, ghost_type); const UInt polynomial_degree = IntegrationOrderFunctor::template getOrder(); quads = GaussIntegrationElement::getQuadraturePoints(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss:: computeJacobianOnQuadPointsByElement(const Matrix & node_coords, const Matrix & quad, Vector & jacobians) const { // jacobian ElementClass::computeJacobian(quad, node_coords, jacobians); } /* -------------------------------------------------------------------------- */ template IntegratorGauss::IntegratorGauss( const Mesh & mesh, const ID & id, const MemoryID & memory_id) : Integrator(mesh, id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void IntegratorGauss::checkJacobians( const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = this->quadrature_points(type, ghost_type).cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); Real * jacobians_val = jacobians(type, ghost_type).storage(); for (UInt i = 0; i < nb_element * nb_quadrature_points; ++i, ++jacobians_val) { if (*jacobians_val < 0) AKANTU_DEBUG_ERROR( "Negative jacobian computed," << " possible problem in the element node ordering (Quadrature Point " << i % nb_quadrature_points << ":" << i / nb_quadrature_points << ":" << type << ":" << ghost_type << ")"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void IntegratorGauss:: computeJacobiansOnIntegrationPoints( const Array & nodes, const Matrix & quad_points, Array & jacobians, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = quad_points.cols(); UInt nb_element = mesh.getNbElement(type, ghost_type); jacobians.resize(nb_element * nb_quadrature_points); auto jacobians_it = jacobians.begin_reinterpret(nb_quadrature_points, nb_element); auto jacobians_begin = jacobians_it; Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type, filter_elements); auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); nb_element = x_el.size(); // Matrix local_coord(spatial_dimension, nb_nodes_per_element); for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) { const Matrix & x = *x_it; Vector & J = *jacobians_it; computeJacobianOnQuadPointsByElement(x, quad_points, J); if (filter_elements == empty_filter) { ++jacobians_it; } else { jacobians_it = jacobians_begin + filter_elements(elem); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) template <> template void IntegratorGauss<_ek_cohesive, DefaultIntegrationOrderFunctor>:: computeJacobiansOnIntegrationPoints( const Array & nodes, const Matrix & quad_points, Array & jacobians, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = quad_points.cols(); UInt nb_element = mesh.getNbElement(type, ghost_type); jacobians.resize(nb_element * nb_quadrature_points); auto jacobians_it = jacobians.begin_reinterpret(nb_quadrature_points, nb_element); auto jacobians_begin = jacobians_it; Vector weights = GaussIntegrationElement::getWeights(); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type, filter_elements); auto x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); UInt nb_nodes_per_subelement = nb_nodes_per_element / 2; Matrix x(spatial_dimension, nb_nodes_per_subelement); nb_element = x_el.size(); - for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) { + auto compute = [&](const auto & el) { + Vector J(jacobians_begin[el]); for (UInt s = 0; s < spatial_dimension; ++s) for (UInt n = 0; n < nb_nodes_per_subelement; ++n) x(s, n) = ((*x_it)(s, n) + (*x_it)(s, n + nb_nodes_per_subelement)) * .5; - Vector & J = *jacobians_it; - if (type == _cohesive_1d_2) J(0) = 1; else - computeJacobianOnQuadPointsByElement(x, quad_points, J); + this->computeJacobianOnQuadPointsByElement(x, quad_points, J); + }; + + for_each_elements(nb_element, filter_elements, compute); - if (filter_elements == empty_filter) { - ++jacobians_it; - } else { - jacobians_it = jacobians_begin + filter_elements(elem); - } - } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ template template void IntegratorGauss:: precomputeJacobiansOnQuadraturePoints(const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Array & jacobians_tmp = jacobians.alloc(0, 1, type, ghost_type); this->computeJacobiansOnIntegrationPoints( nodes, quadrature_points(type, ghost_type), jacobians_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void IntegratorGauss::multiplyJacobiansByWeights( Array & jacobians) const { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_element = jacobians.size() / nb_quadrature_points; Vector weights = GaussIntegrationElement::getWeights(); Array::vector_iterator jacobians_it = jacobians.begin_reinterpret(nb_quadrature_points, nb_element); Array::vector_iterator jacobians_end = jacobians.end_reinterpret(nb_quadrature_points, nb_element); for (; jacobians_it != jacobians_end; ++jacobians_it) { Vector & J = *jacobians_it; J *= weights; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void IntegratorGauss::integrate( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const Array & jacobians, UInt nb_element) const { AKANTU_DEBUG_IN(); intf.resize(nb_element); if (nb_element == 0) return; UInt nb_points = jacobians.size() / nb_element; Array::const_matrix_iterator J_it; Array::matrix_iterator inte_it; Array::const_matrix_iterator f_it; f_it = in_f.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element); inte_it = intf.begin_reinterpret(nb_degree_of_freedom, 1, nb_element); J_it = jacobians.begin_reinterpret(nb_points, 1, nb_element); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Matrix & f = *f_it; const Matrix & J = *J_it; Matrix & inte_f = *inte_it; inte_f.mul(f, J); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void IntegratorGauss::integrate( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); const Array & jac_loc = jacobians(type, ghost_type); if (filter_elements != empty_filter) { UInt nb_element = filter_elements.size(); Array * filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); this->integrate(in_f, intf, nb_degree_of_freedom, *filtered_J, nb_element); delete filtered_J; } else { UInt nb_element = mesh.getNbElement(type, ghost_type); this->integrate(in_f, intf, nb_degree_of_freedom, jac_loc, nb_element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void IntegratorGauss::integrate( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); Matrix quads = this->getIntegrationPoints(); Array jacobians; this->computeJacobiansOnIntegrationPoints(mesh.getNodes(), quads, jacobians, ghost_type); this->multiplyJacobiansByWeights(jacobians); this->integrate(in_f, intf, nb_degree_of_freedom, jacobians, mesh.getNbElement(type, ghost_type)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template Real IntegratorGauss::integrate( const Array & in_f, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); Array intfv(0, 1); integrate(in_f, intfv, 1, ghost_type); Real res = Math::reduce(intfv); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ template template Real IntegratorGauss::integrate( const Array & in_f, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); Array intfv(0, 1); integrate(in_f, intfv, 1, ghost_type, filter_elements); Real res = Math::reduce(intfv); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ template void IntegratorGauss:: integrateOnIntegrationPoints(const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const Array & jacobians, UInt nb_element) const { AKANTU_DEBUG_IN(); UInt nb_points = jacobians.size() / nb_element; Array::const_scalar_iterator J_it; Array::vector_iterator inte_it; Array::const_vector_iterator f_it; intf.resize(nb_element * nb_points); J_it = jacobians.begin(); f_it = in_f.begin(nb_degree_of_freedom); inte_it = intf.begin(nb_degree_of_freedom); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Real & J = *J_it; const Vector & f = *f_it; Vector & inte_f = *inte_it; inte_f = f; inte_f *= J; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void IntegratorGauss:: integrateOnIntegrationPoints(const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); const Array & jac_loc = this->jacobians(type, ghost_type); if (filter_elements != empty_filter) { UInt nb_element = filter_elements.size(); Array * filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); this->integrateOnIntegrationPoints(in_f, intf, nb_degree_of_freedom, *filtered_J, nb_element); } else { UInt nb_element = mesh.getNbElement(type, ghost_type); this->integrateOnIntegrationPoints(in_f, intf, nb_degree_of_freedom, jac_loc, nb_element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void IntegratorGauss::onElementsAdded( - const Array & elements) { + const Array & new_elements) { const auto & nodes = mesh.getNodes(); - for (auto elements_range : MeshElementsByTypes(elements)) { + for (auto elements_range : MeshElementsByTypes(new_elements)) { auto type = elements_range.getType(); auto ghost_type = elements_range.getGhostType(); if (mesh.getKind(type) != kind) continue; auto & elements = elements_range.getElements(); if (not jacobians.exists(type, ghost_type)) jacobians.alloc(0, 1, type, ghost_type); this->computeJacobiansOnIntegrationPoints( nodes, quadrature_points(type, ghost_type), jacobians(type, ghost_type), type, ghost_type, elements); } } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss::initIntegrator( const Array & nodes, const GhostType & ghost_type) { computeQuadraturePoints(ghost_type); precomputeJacobiansOnQuadraturePoints(nodes, ghost_type); checkJacobians(ghost_type); constexpr UInt polynomial_degree = IntegrationOrderFunctor::template getOrder(); multiplyJacobiansByWeights( this->jacobians(type, ghost_type)); } namespace integrator { namespace details { template struct GaussIntegratorInitHelper {}; #define INIT_INTEGRATOR(type) \ _int.template initIntegrator(nodes, ghost_type) #define AKANTU_GAUSS_INTERGRATOR_INIT_HELPER(kind) \ template <> struct GaussIntegratorInitHelper { \ template \ static void call(IntegratorGauss & _int, \ const Array & nodes, const ElementType & type, \ const GhostType & ghost_type) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INIT_INTEGRATOR, kind); \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_GAUSS_INTERGRATOR_INIT_HELPER) #undef AKANTU_GAUSS_INTERGRATOR_INIT_HELPER #undef INIT_INTEGRATOR } // namespace details } // namespace integrator template inline void IntegratorGauss::initIntegrator( const Array & nodes, const ElementType & type, const GhostType & ghost_type) { integrator::details::GaussIntegratorInitHelper::call(*this, nodes, type, ghost_type); } namespace integrator { namespace details { template struct GaussIntegratorComputeJacobiansHelper {}; #define AKANTU_COMPUTE_JACOBIANS(type) \ _int.template computeJacobiansOnIntegrationPoints( \ nodes, quad_points, jacobians, ghost_type, filter_elements); #define AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS(kind) \ template <> struct GaussIntegratorComputeJacobiansHelper { \ template \ static void \ call(const IntegratorGauss & _int, const Array & nodes, \ const Matrix & quad_points, Array & jacobians, \ const ElementType & type, const GhostType & ghost_type, \ const Array & filter_elements) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(AKANTU_COMPUTE_JACOBIANS, kind); \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS) #undef AKANTU_GAUSS_INTERGRATOR_COMPUTE_JACOBIANS #undef AKANTU_COMPUTE_JACOBIANS } // namespace details } // namespace integrator template void IntegratorGauss:: computeJacobiansOnIntegrationPoints( const Array & nodes, const Matrix & quad_points, Array & jacobians, const ElementType & type, const GhostType & ghost_type, const Array & filter_elements) const { integrator::details::GaussIntegratorComputeJacobiansHelper::call( *this, nodes, quad_points, jacobians, type, ghost_type, filter_elements); } } // namespace akantu diff --git a/src/fe_engine/shape_cohesive_inline_impl.cc b/src/fe_engine/shape_cohesive_inline_impl.cc index fd10aa7f3..1717da77c 100644 --- a/src/fe_engine/shape_cohesive_inline_impl.cc +++ b/src/fe_engine/shape_cohesive_inline_impl.cc @@ -1,310 +1,301 @@ /** * @file shape_cohesive_inline_impl.cc * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Feb 03 2012 * @date last modification: Thu Oct 15 2015 * * @brief ShapeCohesive 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 "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_COHESIVE_INLINE_IMPL_CC__ #define __AKANTU_SHAPE_COHESIVE_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline ShapeLagrange<_ek_cohesive>::ShapeLagrange(const Mesh & mesh, const ID & id, const MemoryID & memory_id) : ShapeLagrangeBase(mesh, _ek_cohesive, id, memory_id) {} #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); /* -------------------------------------------------------------------------- */ inline void ShapeLagrange<_ek_cohesive>::initShapeFunctions( const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type) { AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeShapeDerivativesOnIntegrationPoints( const Array &, 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 size_of_shapesd = ElementClass::getShapeDerivativesSize(); + UInt spatial_dimension = ElementClass::getNaturalSpaceDimension(); + 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); - Real * shapesd_val = shape_derivatives.storage(); - if (filter_elements == empty_filter) - nb_element = filter_elements.size(); - - for (UInt elem = 0; elem < nb_element; ++elem) { - Tensor3 B(shapesd_val, spatial_dimension, nb_nodes_per_element, - nb_points); + auto compute = [&](auto & el) { + auto ptr = shapesd_val + el * nb_points * size_of_shapesd; + Tensor3 B(ptr, spatial_dimension, nb_nodes_per_element, nb_points); ElementClass::computeDNDS(integration_points, B); + }; - if (filter_elements != empty_filter) - shapesd_val += size_of_shapesd * nb_points; - else { - shapesd_val = shape_derivatives.storage() + - filter_elements(elem) * size_of_shapesd * nb_points; - } - } + for_each_elements(nb_element, filter_elements, compute); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void ShapeLagrange<_ek_cohesive>::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_COHESIVE_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES); #undef AKANTU_COMPUTE_SHAPES } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::precomputeShapesOnIntegrationPoints( const Array & nodes, 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 void ShapeLagrange<_ek_cohesive>::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, 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 void ShapeLagrange<_ek_cohesive>::extractNodalToElementField( const Array & nodal_f, Array & elemental_f, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_itp_element = ElementClass::getNbNodesPerInterpolationElement(); - UInt nb_nodes_per_element = ElementClass::getNbNodesPerElement(); UInt nb_degree_of_freedom = nodal_f.getNbComponent(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); - UInt * conn_val = this->mesh.getConnectivity(type, ghost_type).storage(); + + const auto & conn_array = this->mesh.getConnectivity(type, ghost_type); + auto conn = conn_array.begin(conn_array.getNbComponent() / 2, 2); if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } elemental_f.resize(nb_element); Array::matrix_iterator u_it = elemental_f.begin(nb_degree_of_freedom, nb_nodes_per_itp_element); - UInt * el_conn; ReduceFunction reduce_function; - for (UInt el = 0; el < nb_element; ++el, ++u_it) { + auto compute = [&](auto & el) { Matrix & u = *u_it; - if (filter_elements != empty_filter) - el_conn = conn_val + filter_elements(el) * nb_nodes_per_element; - else - el_conn = conn_val + el * nb_nodes_per_element; + Matrix el_conn(conn[el]); // compute the average/difference of the nodal field loaded from cohesive // element - for (UInt n = 0; n < nb_nodes_per_itp_element; ++n) { - UInt node_plus = *(el_conn + n); - UInt node_minus = *(el_conn + n + nb_nodes_per_itp_element); + for (UInt n = 0; n < el_conn.rows(); ++n) { + UInt node_plus = el_conn(n, 0); + UInt node_minus = el_conn(n, 1); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { Real u_plus = nodal_f(node_plus, d); Real u_minus = nodal_f(node_minus, d); u(d, n) = reduce_function(u_plus, u_minus); } } - } + + ++u_it; + }; + + for_each_elements(nb_element, filter_elements, compute); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::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(this->shapes.exists(itp_type, ghost_type), "No shapes for the type " << this->shapes.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); this->extractNodalToElementField(in_u, u_el, ghost_type, filter_elements); this->template interpolateElementalFieldOnIntegrationPoints( u_el, out_uq, ghost_type, shapes(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::variationOnIntegrationPoints( const Array & in_u, Array & 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( this->shapes_derivatives.exists(itp_type, ghost_type), "No shapes for the type " << this->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); this->extractNodalToElementField(in_u, u_el, ghost_type, filter_elements); this->template gradientElementalFieldOnIntegrationPoints( u_el, nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeNormalsOnIntegrationPoints( const Array & u, Array & normals_u, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); UInt nb_points = this->integration_points(type, ghost_type).cols(); UInt spatial_dimension = this->mesh.getSpatialDimension(); if (filter_elements != empty_filter) nb_element = filter_elements.size(); normals_u.resize(nb_points * nb_element); Array tangents_u(nb_element * nb_points, (spatial_dimension * (spatial_dimension - 1))); this->template variationOnIntegrationPoints( u, tangents_u, spatial_dimension, ghost_type, filter_elements); Array::vector_iterator normal = normals_u.begin(spatial_dimension); Array::vector_iterator normal_end = normals_u.end(spatial_dimension); Real * tangent = tangents_u.storage(); if (spatial_dimension == 3) for (; normal != normal_end; ++normal) { Math::vectorProduct3(tangent, tangent + spatial_dimension, normal->storage()); (*normal) /= normal->norm(); tangent += spatial_dimension * 2; } else if (spatial_dimension == 2) for (; normal != normal_end; ++normal) { Vector a1(tangent, spatial_dimension); (*normal)(0) = -a1(1); (*normal)(1) = a1(0); (*normal) /= normal->norm(); tangent += spatial_dimension; } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* __AKANTU_SHAPE_COHESIVE_INLINE_IMPL_CC__ */ diff --git a/src/fe_engine/shape_lagrange_base_inline_impl.cc b/src/fe_engine/shape_lagrange_base_inline_impl.cc index e9dd7163b..1c09adcef 100644 --- a/src/fe_engine/shape_lagrange_base_inline_impl.cc +++ b/src/fe_engine/shape_lagrange_base_inline_impl.cc @@ -1,97 +1,97 @@ /** * @file shape_lagrange_base_inline_impl.cc * * @author Nicolas Richart * * @date creation Thu Jul 27 2017 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "shape_lagrange_base.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_LAGRANGE_BASE_INLINE_IMPL_CC__ #define __AKANTU_SHAPE_LAGRANGE_BASE_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline const Array & ShapeLagrangeBase::getShapes(const ElementType & el_type, const GhostType & ghost_type) const { return shapes(FEEngine::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ inline const Array & ShapeLagrangeBase::getShapesDerivatives(const ElementType & el_type, const GhostType & ghost_type) const { return shapes_derivatives(FEEngine::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ template void ShapeLagrangeBase::computeShapesOnIntegrationPoints( const Array &, const Matrix & integration_points, Array & shapes, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); shapes.resize(nb_element * nb_points); #if !defined(AKANTU_NDEBUG) UInt size_of_shapes = ElementClass::getShapeSize(); 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( ElementClass::getNbNodesPerInterpolationElement(), nb_points, nb_element); auto shapes_begin = shapes_it; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } for (UInt elem = 0; elem < nb_element; ++elem) { + if (filter_elements != empty_filter) + shapes_it = shapes_begin + filter_elements(elem); + Matrix & N = *shapes_it; ElementClass::computeShapes(integration_points, N); - if (filter_elements == empty_filter) { + if (filter_elements == empty_filter) ++shapes_it; - } else { - shapes_it = shapes_begin + filter_elements(elem); - } } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* __AKANTU_SHAPE_LAGRANGE_BASE_INLINE_IMPL_CC__ */ diff --git a/src/fe_engine/shape_lagrange_inline_impl.cc b/src/fe_engine/shape_lagrange_inline_impl.cc index 367a894fd..cd243ae77 100644 --- a/src/fe_engine/shape_lagrange_inline_impl.cc +++ b/src/fe_engine/shape_lagrange_inline_impl.cc @@ -1,411 +1,411 @@ /** * @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 "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 dxds Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); // compute shape derivatives 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; - else { - shapesd_val = shape_derivatives.storage() + - filter_elements(elem) * 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::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->mul(*field_it, *shapes_it); } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* __AKANTU_SHAPE_LAGRANGE_INLINE_IMPL_CC__ */ diff --git a/src/mesh_utils/cohesive_element_inserter.cc b/src/mesh_utils/cohesive_element_inserter.cc index af9656d93..d22c3f001 100644 --- a/src/mesh_utils/cohesive_element_inserter.cc +++ b/src/mesh_utils/cohesive_element_inserter.cc @@ -1,458 +1,448 @@ /** * @file cohesive_element_inserter.cc * * @author Marco Vocialta * * @date creation: Wed Dec 04 2013 * @date last modification: Sun Oct 04 2015 * * @brief Cohesive element inserter functions * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "cohesive_element_inserter.hh" #include "element_group.hh" #include #include /* -------------------------------------------------------------------------- */ namespace akantu { CohesiveElementInserter::CohesiveElementInserter(Mesh & mesh, bool is_extrinsic, const ID & id) : id(id), mesh(mesh), mesh_facets(mesh.initMeshFacets()), insertion_facets("insertion_facets", id), insertion_limits(mesh.getSpatialDimension(), 2), check_facets("check_facets", id) { MeshUtils::buildAllFacets(mesh, mesh_facets, 0); init(is_extrinsic); } /* -------------------------------------------------------------------------- */ CohesiveElementInserter::~CohesiveElementInserter() { #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) delete global_ids_updater; #endif } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::init(bool is_extrinsic) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); MeshUtils::resetFacetToDouble(mesh_facets); /// initialize facet insertion array insertion_facets.initialize(mesh_facets, _spatial_dimension = (spatial_dimension - 1), _with_nb_element = true); // mesh_facets.initElementTypeMapArray( // insertion_facets, 1, spatial_dimension - 1, false, _ek_regular, true); /// init insertion limits for (UInt dim = 0; dim < spatial_dimension; ++dim) { insertion_limits(dim, 0) = std::numeric_limits::max() * (-1.); insertion_limits(dim, 1) = std::numeric_limits::max(); } if (is_extrinsic) { check_facets.initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(check_facets, 1, spatial_dimension - // 1); initFacetsCheck(); } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) facet_synchronizer = NULL; global_ids_updater = NULL; #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::initFacetsCheck() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType facet_gt = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, facet_gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, facet_gt); for (; it != last; ++it) { ElementType facet_type = *it; Array & f_check = check_facets(facet_type, facet_gt); const Array> & element_to_facet = mesh_facets.getElementToSubelement(facet_type, facet_gt); UInt nb_facet = element_to_facet.size(); f_check.resize(nb_facet); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull || (element_to_facet(f)[0].ghost_type == _ghost && element_to_facet(f)[1].ghost_type == _ghost)) { f_check(f) = false; } else f_check(f) = true; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::limitCheckFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector bary_facet(spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); for (; it != end; ++it) { ElementType type = *it; Array & f_check = check_facets(type, ghost_type); UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); for (UInt f = 0; f < nb_facet; ++f) { if (f_check(f)) { mesh_facets.getBarycenter(f, type, bary_facet.storage(), ghost_type); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit != spatial_dimension) f_check(f) = false; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::setLimit(SpacialDirection axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_ASSERT( axis < mesh.getSpatialDimension(), "You are trying to limit insertion in a direction that doesn't exist"); insertion_limits(axis, 0) = std::min(first_limit, second_limit); insertion_limits(axis, 1) = std::max(first_limit, second_limit); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertIntrinsicElements() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); Vector bary_facet(spatial_dimension); - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - - GhostType ghost_type = *gt; - - Mesh::type_iterator it = - mesh_facets.firstType(spatial_dimension - 1, ghost_type); - Mesh::type_iterator end = - mesh_facets.lastType(spatial_dimension - 1, ghost_type); - - for (; it != end; ++it) { - const ElementType type_facet = *it; - Array & f_insertion = insertion_facets(type_facet, ghost_type); - Array> & element_to_facet = + for (auto && ghost_type : ghost_types) { + for (const auto & type_facet : + mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { + auto & f_insertion = insertion_facets(type_facet, ghost_type); + auto & element_to_facet = mesh_facets.getElementToSubelement(type_facet, ghost_type); UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type); - for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull) continue; mesh_facets.getBarycenter(f, type_facet, bary_facet.storage(), ghost_type); UInt coord_in_limit = 0; while (coord_in_limit < spatial_dimension && bary_facet(coord_in_limit) > insertion_limits(coord_in_limit, 0) && bary_facet(coord_in_limit) < insertion_limits(coord_in_limit, 1)) ++coord_in_limit; if (coord_in_limit == spatial_dimension) f_insertion(f) = true; } } } AKANTU_DEBUG_OUT(); return insertElements(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::insertIntrinsicElements(std::string physname, UInt material_index) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray * phys_data; try { phys_data = &(mesh_facets.getData("physical_names")); } catch (...) { phys_data = &(mesh_facets.registerData("physical_names")); phys_data->initialize(mesh_facets, _spatial_dimension = spatial_dimension - 1, _with_nb_element = true); // mesh_facets.initElementTypeMapArray(*phys_data, 1, spatial_dimension - 1, // false, _ek_regular, true); } Vector bary_facet(spatial_dimension); mesh_facets.createElementGroup(physname); GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); for (; it != end; ++it) { const ElementType type_facet = *it; Array & f_insertion = insertion_facets(type_facet, ghost_type); Array> & element_to_facet = mesh_facets.getElementToSubelement(type_facet, ghost_type); UInt nb_facet = mesh_facets.getNbElement(type_facet, ghost_type); UInt coord_in_limit = 0; ElementGroup & group = mesh.getElementGroup(physname); ElementGroup & group_facet = mesh_facets.getElementGroup(physname); Vector bary_physgroup(spatial_dimension); Real norm_bary; for (ElementGroup::const_element_iterator el_it( group.begin(type_facet, ghost_type)); el_it != group.end(type_facet, ghost_type); ++el_it) { UInt e = *el_it; mesh.getBarycenter(e, type_facet, bary_physgroup.storage(), ghost_type); #ifndef AKANTU_NDEBUG bool find_a_partner = false; #endif norm_bary = bary_physgroup.norm(); Array & material_id = (*phys_data)(type_facet, ghost_type); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)[1] == ElementNull) continue; mesh_facets.getBarycenter(f, type_facet, bary_facet.storage(), ghost_type); coord_in_limit = 0; while (coord_in_limit < spatial_dimension && (std::abs(bary_facet(coord_in_limit) - bary_physgroup(coord_in_limit)) / norm_bary < Math::getTolerance())) ++coord_in_limit; if (coord_in_limit == spatial_dimension) { f_insertion(f) = true; #ifndef AKANTU_NDEBUG find_a_partner = true; #endif group_facet.add(type_facet, f, ghost_type, false); material_id(f) = material_index; break; } } AKANTU_DEBUG_ASSERT(find_a_partner, "The element nO " << e << " of physical group " << physname << " did not find its associated facet!" << " Try to decrease math tolerance. " << std::endl); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt CohesiveElementInserter::insertElements(bool only_double_facets) { CohesiveNewNodesEvent node_event; NewElementsEvent element_event; Array new_pairs(0, 2); UInt nb_new_elements = MeshUtils::insertCohesiveElements( - mesh, mesh_facets, insertion_facets, new_pairs, - element_event.getList(), only_double_facets); + mesh, mesh_facets, insertion_facets, new_pairs, element_event.getList(), + only_double_facets); UInt nb_new_nodes = new_pairs.size(); node_event.getList().reserve(nb_new_nodes); node_event.getOldNodesList().reserve(nb_new_nodes); - for(UInt n = 0; n < nb_new_nodes; ++n) { + for (UInt n = 0; n < nb_new_nodes; ++n) { node_event.getList().push_back(new_pairs(n, 1)); node_event.getOldNodesList().push_back(new_pairs(n, 0)); } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (mesh.getNodesType().size()) { /// update nodes type updateNodesType(mesh, node_event); updateNodesType(mesh_facets, node_event); /// update global ids nb_new_nodes = this->updateGlobalIDs(node_event); /// compute total number of new elements StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); comm.allReduce(&nb_new_elements, 1, _so_sum); } #endif if (nb_new_nodes > 0) mesh.sendEvent(node_event); if (nb_new_elements > 0) { updateInsertionFacets(); - //mesh.updateTypesOffsets(_not_ghost); + // mesh.updateTypesOffsets(_not_ghost); mesh.sendEvent(element_event); MeshUtils::resetFacetToDouble(mesh_facets); } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (mesh.getNodesType().size()) { /// update global ids this->synchronizeGlobalIDs(node_event); } #endif return nb_new_elements; } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::updateInsertionFacets() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType facet_gt = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, facet_gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, facet_gt); for (; it != last; ++it) { ElementType facet_type = *it; Array & ins_facets = insertion_facets(facet_type, facet_gt); // this is the extrinsic case if (check_facets.exists(facet_type, facet_gt)) { Array & f_check = check_facets(facet_type, facet_gt); UInt nb_facets = f_check.size(); for (UInt f = 0; f < ins_facets.size(); ++f) { if (ins_facets(f)) { ++nb_facets; ins_facets(f) = false; f_check(f) = false; } } f_check.resize(nb_facets); } // and this the intrinsic one else { ins_facets.resize(mesh_facets.getNbElement(facet_type, facet_gt)); ins_facets.set(false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CohesiveElementInserter::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "CohesiveElementInserter [" << std::endl; stream << space << " + mesh [" << std::endl; mesh.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + mesh_facets [" << std::endl; mesh_facets.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index 06e8d6139..5ebab6973 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,788 +1,775 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Jan 13 2016 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "dumpable_inline_impl.hh" #include "material_cohesive.hh" #include "shape_cohesive.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); inserter = NULL; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) facet_synchronizer = NULL; facet_stress_synchronizer = NULL; cohesive_element_synchronizer = NULL; global_connectivity = NULL; #endif delete material_selector; material_selector = new DefaultMaterialCohesiveSelector(*this); this->registerEventHandler(*this, _ehp_solid_mechanics_model_cohesive); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() { AKANTU_DEBUG_IN(); delete inserter; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) delete cohesive_element_synchronizer; delete facet_synchronizer; delete facet_stress_synchronizer; #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step) { SolidMechanicsModel::setTimeStep(time_step); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFull(const ModelOptions & options) { AKANTU_DEBUG_IN(); const SolidMechanicsModelCohesiveOptions & smmc_options = dynamic_cast(options); this->is_extrinsic = smmc_options.extrinsic; if (!inserter) inserter = new CohesiveElementInserter(mesh, is_extrinsic, id + ":cohesive_element_inserter"); SolidMechanicsModel::initFull(options); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if (!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt cohesive_index = 0; while ((dynamic_cast(materials[cohesive_index].get()) == nullptr) && cohesive_index <= materials.size()) ++cohesive_index; AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(), "No cohesive materials in the material input file"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion if (is_extrinsic) { const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize(mesh_facets, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(facet_material, 1, // spatial_dimension - 1); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : mesh_facets.elementTypes(Model::spatial_dimension - 1, ghost_type)) { element.type = type; Array & f_material = facet_material(type, ghost_type); UInt nb_element = mesh_facets.getNbElement(type, ghost_type); f_material.resize(nb_element); f_material.set(cohesive_index); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt mat_index = (*material_selector)(element); f_material(el) = mat_index; MaterialCohesive & mat = dynamic_cast(*materials[mat_index]); mat.addFacet(element); } } } SolidMechanicsModel::initMaterials(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif initAutomaticInsertion(); } else { // TODO think of something a bit mor consistant than just coding the first // thing that comes in Fabian's head.... typedef ParserSection::const_section_iterator const_section_iterator; std::pair sub_sections = this->parser->getSubSections(_st_mesh); if (sub_sections.first != sub_sections.second) { std::string cohesive_surfaces = sub_sections.first->getParameter("cohesive_surfaces"); this->initIntrinsicCohesiveMaterials(cohesive_surfaces); } else { this->initIntrinsicCohesiveMaterials(cohesive_index); } } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( std::string cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif std::istringstream split(cohesive_surfaces); std::string physname; while (std::getline(split, physname, ',')) { AKANTU_DEBUG_INFO( "Pre-inserting cohesive elements along facets from physical surface: " << physname); insertElementsFromMeshData(physname); } synchronizeInsertionData(); SolidMechanicsModel::initMaterials(); if (is_default_material_selector) delete material_selector; material_selector = new MeshDataMaterialCohesiveSelector(*this); // UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeInsertionData() { #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) { facet_synchronizer->asynchronousSynchronize(*inserter, _gst_ce_groups); facet_synchronizer->waitEndSynchronize(*inserter, _gst_ce_groups); } #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( UInt cohesive_index) { AKANTU_DEBUG_IN(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(Model::spatial_dimension, *gt, _ek_cohesive); Mesh::type_iterator last = mesh.lastType(Model::spatial_dimension, *gt, _ek_cohesive); for (; first != last; ++first) { Array & mat_indexes = this->material_index(*first, *gt); Array & mat_loc_num = this->material_local_numbering(*first, *gt); mat_indexes.set(cohesive_index); mat_loc_num.clear(); } } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("CohesiveFEEngine", mesh, Model::spatial_dimension); /// add cohesive type connectivity ElementType type = _not_defined; - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - - GhostType type_ghost = *gt; + for (auto && type_ghost : ghost_types) { + for (const auto & tmp_type : mesh.elementTypes(spatial_dimension, type_ghost)) { + const Array & connectivity = mesh.getConnectivity(tmp_type, type_ghost); + if (connectivity.size() == 0) + continue; - Mesh::type_iterator it = - mesh.firstType(Model::spatial_dimension, type_ghost); - Mesh::type_iterator last = - mesh.lastType(Model::spatial_dimension, type_ghost); - - for (; it != last; ++it) { - const Array & connectivity = mesh.getConnectivity(*it, type_ghost); - if (connectivity.size() != 0) { - type = *it; - ElementType type_facet = Mesh::getFacetType(type); - ElementType type_cohesive = - FEEngine::getCohesiveElementType(type_facet); - mesh.addConnectivityType(type_cohesive, type_ghost); - } + type = tmp_type; + ElementType type_facet = Mesh::getFacetType(type); + ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); + mesh.addConnectivityType(type_cohesive, type_ghost); } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); registerFEEngineObject( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::limitInsertion(BC::Axis axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_IN(); inserter->setLimit(axis, first_limit, second_limit); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); // UInt nb_new_elements = inserter->insertIntrinsicElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertElementsFromMeshData( std::string physname) { AKANTU_DEBUG_IN(); UInt material_index = SolidMechanicsModel::getMaterialIndex(physname); inserter->insertIntrinsicElements(physname, material_index); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer( *inserter, rank_to_element, *data_accessor); } #endif facet_stress.initialize(inserter->getMeshFacets(), _nb_component = 2 * Model::spatial_dimension * Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // inserter->getMeshFacets().initElementTypeMapArray( // facet_stress, 2 * spatial_dimension * spatial_dimension, // spatial_dimension - 1); resizeFacetStress(); /// compute normals on facets computeNormals(); initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); inserter->limitCheckFacets(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer( *inserter, rank_to_element, *data_accessor); } #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element Array & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); Array & el_q_facet = elements_quad_facets(type, elem_gt); ElementType facet_type = Mesh::getFacetType(type); UInt nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(facet_type); el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element global_facet_elem = facet_to_element(el, f); UInt global_facet = global_facet_elem.element; GhostType facet_gt = global_facet_elem.ghost_type; const Array & quad_f = quad_facets(facet_type, facet_gt); for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < Model::spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_f(global_facet * nb_quad_per_facet + q, s); } } } } } } /// loop over non cohesive materials for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*materials[m]); } catch (std::bad_cast &) { /// initialize the interpolation function materials[m]->initElementalFieldInterpolation(elements_quad_facets); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { MaterialCohesive & mat = dynamic_cast(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); if (isDefaultSolverExplicit()) { for (auto & material : materials) { try { MaterialCohesive & mat = dynamic_cast(*material); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*material); } catch (std::bad_cast &) { /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } } #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData( debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif this->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { interpolateStress(); for (auto & mat : materials) { try { MaterialCohesive & mat_cohesive = dynamic_cast(*mat); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch (std::bad_cast &) { } } /// communicate data among processors this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) updateCohesiveSynchronizers(); #endif SolidMechanicsModel::onElementsAdded(element_list, event); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (cohesive_element_synchronizer != NULL) cohesive_element_synchronizer->computeAllBufferSizes(*this); #endif if (is_extrinsic) resizeFacetStress(); /// if (method != _explicit_lumped_mass) { /// this->initSolver(); /// } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); // Array nodes_list(nb_new_nodes); // for (UInt n = 0; n < nb_new_nodes; ++n) // nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(new_nodes, event); UInt new_node, old_node; try { const auto & cohesive_event = dynamic_cast(event); const auto & old_nodes = cohesive_event.getOldNodesList(); auto copy = [this, &new_node, &old_node](auto & arr) { for (UInt s = 0; s < spatial_dimension; ++s) { arr(new_node, s) = arr(old_node, s); } }; for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; copy(*displacement); - copy(*velocity); - copy(*acceleration); copy(*blocked_dofs); + if (velocity) + copy(*velocity); + + if (acceleration) + copy(*acceleration); + if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); } // if (this->getDOFManager().hasMatrix("M")) { // this->assembleMass(old_nodes); // } // if (this->getDOFManager().hasLumpedMatrix("M")) { // this->assembleMassLumped(old_nodes); // } } catch (std::bad_cast &) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod &) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ for (auto & mat : materials) { if (mat->isFiniteDeformation()) mat->computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = inserter->getMeshFacets(); - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - GhostType ghost_type = *gt; - - Mesh::type_iterator it = - mesh_facets.firstType(Model::spatial_dimension - 1, ghost_type); - Mesh::type_iterator end = - mesh_facets.lastType(Model::spatial_dimension - 1, ghost_type); - for (; it != end; ++it) { - ElementType type = *it; - + for (auto && ghost_type : ghost_types) { + for (const auto & type : + mesh_facets.elementTypes(spatial_dimension - 1, ghost_type)) { UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); UInt nb_quadrature_points = getFEEngine("FacetsFEEngine") .getNbIntegrationPoints(type, ghost_type); UInt new_size = nb_facet * nb_quadrature_points; facet_stress(type, ghost_type).resize(new_size); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh index dbd59372a..54952cde6 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh @@ -1,297 +1,300 @@ /** * @file solid_mechanics_model_cohesive.hh * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Mon Dec 14 2015 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ #include "cohesive_element_inserter.hh" #include "material_selector_cohesive.hh" #include "solid_mechanics_model.hh" #include "solid_mechanics_model_event_handler.hh" #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) #include "facet_stress_synchronizer.hh" #include "facet_synchronizer.hh" #endif + +#include "random_internal_field.hh" // included to have the specialization of + // ParameterTyped::operator Real() /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ struct FacetsCohesiveIntegrationOrderFunctor { template ::cohesive_type> struct _helper { static constexpr int get() { return ElementClassProperty::polynomial_degree; } }; template struct _helper { static constexpr int get() { return ElementClassProperty::polynomial_degree; } }; template static inline constexpr int getOrder() { return _helper::get(); } }; /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelCohesiveOptions : public SolidMechanicsModelOptions { SolidMechanicsModelCohesiveOptions( AnalysisMethod analysis_method = _explicit_lumped_mass, bool extrinsic = false) : SolidMechanicsModelOptions(analysis_method), extrinsic(extrinsic) {} template SolidMechanicsModelCohesiveOptions(use_named_args_t, pack &&... _pack) - : SolidMechanicsModelOptions( + : SolidMechanicsModelCohesiveOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass), OPTIONAL_NAMED_ARG(is_extrinsic, false)) {} bool extrinsic{false}; }; /* -------------------------------------------------------------------------- */ /* Solid Mechanics Model for Cohesive elements */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModelCohesive : public SolidMechanicsModel, public SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewCohesiveNodesEvent : public NewNodesEvent { public: AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array &); AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array &); protected: Array old_nodes; }; typedef FEEngineTemplate MyFEEngineCohesiveType; typedef FEEngineTemplate MyFEEngineFacetType; SolidMechanicsModelCohesive(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model_cohesive", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModelCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step); /// assemble the residual for the explicit scheme virtual void assembleInternalForces(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function to perform a stress check on each facet and insert /// cohesive elements if needed (returns the number of new cohesive /// elements) UInt checkCohesiveStress(); /// interpolate stress on facets void interpolateStress(); /// initialize the cohesive model void initFull(const ModelOptions & options = SolidMechanicsModelCohesiveOptions()); template void initFull(named_argument::param_t && first, pack &&... _pack) { this->initFull( SolidMechanicsModelCohesiveOptions{use_named_args, first, _pack...}); } /// initialize the model void initModel(); /// initialize cohesive material void initMaterials(); /// init facet filters for cohesive materials void initFacetFilter(); /// limit the cohesive element insertion to a given area void limitInsertion(BC::Axis axis, Real first_limit, Real second_limit); /// update automatic insertion after a change in the element inserter void updateAutomaticInsertion(); /// insert intrinsic cohesive elements void insertIntrinsicElements(); // synchronize facets physical data before insertion along physical surfaces void synchronizeInsertionData(); // template bool solveStepCohesive(Real tolerance, Real & error, UInt // max_iteration = 100, // bool load_reduction = false, // Real tol_increase_factor = 1.0, // bool do_not_factorize = false); /// initialize stress interpolation void initStressInterpolation(); private: /// initialize cohesive material with intrinsic insertion (by default) void initIntrinsicCohesiveMaterials(UInt cohesive_index); /// initialize cohesive material with intrinsic insertion (if physical /// surfaces are precised) void initIntrinsicCohesiveMaterials(std::string cohesive_surfaces); /// insert cohesive elements along a given physical surface of the mesh void insertElementsFromMeshData(std::string physical_name); /// initialize completely the model for extrinsic elements void initAutomaticInsertion(); /// compute facets' normals void computeNormals(); /// resize facet stress void resizeFacetStress(); /// init facets_check array void initFacetsCheck(); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event); virtual void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event); /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void onEndSolveStep(const AnalysisMethod & method); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get facet mesh AKANTU_GET_MACRO(MeshFacets, mesh.getMeshFacets(), const Mesh &); /// get stress on facets vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressOnFacets, facet_stress, Real); /// get facet material AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetMaterial, facet_material, UInt); /// get facet material AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetMaterial, facet_material, UInt); /// get facet material AKANTU_GET_MACRO(FacetMaterial, facet_material, const ElementTypeMapArray &); /// @todo THIS HAS TO BE CHANGED AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Tangents, tangents, Real); /// get element inserter AKANTU_GET_MACRO_NOT_CONST(ElementInserter, *inserter, CohesiveElementInserter &); /// get is_extrinsic boolean AKANTU_GET_MACRO(IsExtrinsic, is_extrinsic, bool); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// @todo store tangents when normals are computed: ElementTypeMapArray tangents; /// stress on facets on the two sides by quadrature point ElementTypeMapArray facet_stress; /// material to use if a cohesive element is created on a facet ElementTypeMapArray facet_material; bool is_extrinsic; /// cohesive element inserter CohesiveElementInserter * inserter; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive_parallel.hh" #endif }; /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const SolidMechanicsModelCohesive & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "solid_mechanics_model_cohesive_inline_impl.cc" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc index 479f7a3a9..ba571c627 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc @@ -1,191 +1,190 @@ /** * @file test_cohesive_buildfragments.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Thu Oct 15 2015 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "fragment_manager.hh" #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); Math::setTolerance(1e-14); const UInt spatial_dimension = 2; const UInt max_steps = 200; Real strain_rate = 1.e5; ElementType type = _quadrangle_4; Real L = 0.03; Real theoretical_mass = L * L / 20. * 2500; ElementType type_facet = Mesh::getFacetType(type); ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); Mesh mesh(spatial_dimension); mesh.read("mesh.msh"); - mesh.createGroupsFromMeshData("physical_names"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull( SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); Real time_step = model.getStableTimeStep() * 0.05; model.setTimeStep(time_step); // std::cout << "Time step: " << time_step << std::endl; Real disp_increment = strain_rate * L / 2. * time_step; model.assembleMassLumped(); Array & velocity = model.getVelocity(); const Array & position = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); /// initial conditions for (UInt n = 0; n < nb_nodes; ++n) velocity(n, 0) = strain_rate * position(n, 0); /// boundary conditions model.applyBC(BC::Dirichlet::FixedValue(0, _x), "Left_side"); model.applyBC(BC::Dirichlet::FixedValue(0, _x), "Right_side"); UInt cohesive_index = 1; UInt nb_quad_per_facet = model.getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); MaterialCohesive & mat_cohesive = dynamic_cast(model.getMaterial(cohesive_index)); const Array & damage = mat_cohesive.getDamage(type_cohesive); FragmentManager fragment_manager(model, false); const Array & fragment_mass = fragment_manager.getMass(); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.checkCohesiveStress(); model.solveStep(); /// apply boundary conditions model.applyBC(BC::Dirichlet::IncrementValue(-disp_increment, _x), "Left_side"); model.applyBC(BC::Dirichlet::IncrementValue(disp_increment, _x), "Right_side"); if (s % 1 == 0) { // model.dump(); std::cout << "passing step " << s << "/" << max_steps << std::endl; fragment_manager.computeAllData(); /// check number of fragments UInt nb_fragment_num = fragment_manager.getNbFragment(); UInt nb_cohesive_elements = mesh.getNbElement(type_cohesive); UInt nb_fragment = 1; for (UInt el = 0; el < nb_cohesive_elements; ++el) { UInt q = 0; while (q < nb_quad_per_facet && Math::are_float_equal(damage(el * nb_quad_per_facet + q), 1)) ++q; if (q == nb_quad_per_facet) { ++nb_fragment; } } if (nb_fragment != nb_fragment_num) { std::cout << "The number of fragments is wrong!" << std::endl; return EXIT_FAILURE; } /// check mass computation Real total_mass = 0.; for (UInt frag = 0; frag < nb_fragment_num; ++frag) { total_mass += fragment_mass(frag); } if (!Math::are_float_equal(theoretical_mass, total_mass)) { std::cout << "The fragments' mass is wrong!" << std::endl; return EXIT_FAILURE; } } } model.dump(); /// check velocities UInt nb_fragment = fragment_manager.getNbFragment(); const Array & fragment_velocity = fragment_manager.getVelocity(); const Array & fragment_center = fragment_manager.getCenterOfMass(); Real fragment_length = L / nb_fragment; Real initial_position = -L / 2. + fragment_length / 2.; for (UInt frag = 0; frag < nb_fragment; ++frag) { Real theoretical_center = initial_position + fragment_length * frag; if (!Math::are_float_equal(fragment_center(frag, 0), theoretical_center)) { std::cout << "The fragments' center is wrong!" << std::endl; return EXIT_FAILURE; } Real initial_vel = fragment_center(frag, 0) * strain_rate; Math::setTolerance(100); if (!Math::are_float_equal(fragment_velocity(frag), initial_vel)) { std::cout << "The fragments' velocity is wrong!" << std::endl; return EXIT_FAILURE; } } finalize(); std::cout << "OK: test_cohesive_buildfragments was passed!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc index 5ab2f4031..b45284d9e 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc @@ -1,138 +1,137 @@ /** * @file test_cohesive_extrinsic.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Mar 18 2015 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); - debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 2; const UInt max_steps = 1000; Mesh mesh(spatial_dimension); mesh.read("triangle.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull( SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); model.limitInsertion(_y, -0.30, -0.20); model.updateAutomaticInsertion(); mesh.setBaseName("test_cohesive_extrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("mass"); model.addDumpField("velocity"); model.addDumpField("acceleration"); - model.addDumpField("external_force"); - model.addDumpField("internal_force"); + model.addDumpFieldVector("external_force"); + model.addDumpFieldVector("internal_force"); model.addDumpField("grad_u"); model.dump(); Real time_step = model.getStableTimeStep() * 0.05; model.setTimeStep(time_step); std::cout << "Time step: " << time_step << std::endl; model.assembleMassLumped(); Array & position = mesh.getNodes(); Array & velocity = model.getVelocity(); Array & boundary = model.getBlockedDOFs(); Array & displacement = model.getDisplacement(); // const Array & residual = model.getResidual(); UInt nb_nodes = mesh.getNbNodes(); /// boundary conditions for (UInt n = 0; n < nb_nodes; ++n) { if (position(n, 1) > 0.99 || position(n, 1) < -0.99) boundary(n, 1) = true; if (position(n, 0) > 0.99 || position(n, 0) < -0.99) boundary(n, 0) = true; } /// initial conditions Real loading_rate = 0.5; Real disp_update = loading_rate * time_step; for (UInt n = 0; n < nb_nodes; ++n) { velocity(n, 1) = loading_rate * position(n, 1); } /// Main loop for (UInt s = 1; s <= max_steps; ++s) { /// update displacement on extreme nodes for (UInt n = 0; n < mesh.getNbNodes(); ++n) { if (position(n, 1) > 0.99 || position(n, 1) < -0.99) displacement(n, 1) += disp_update * position(n, 1); } model.checkCohesiveStress(); model.solveStep(); if (s % 100 == 0) { std::cout << "passing step " << s << "/" << max_steps << std::endl; } model.dump(); } Real Ed = model.getEnergy("dissipated"); Real Edt = 200 * std::sqrt(2); std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_extrinsic was passed!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_tetrahedron.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_tetrahedron.cc index 8186b9595..0b5d70779 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_tetrahedron.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic_tetrahedron.cc @@ -1,238 +1,240 @@ /** * @file test_cohesive_extrinsic_tetrahedron.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Thu Oct 15 2015 * * @brief Test for serial extrinsic cohesive elements for tetrahedron * * @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 #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; Real function(Real constant, Real x, Real y, Real z) { return constant + 2. * x + 3. * y + 4 * z; } int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); // const UInt max_steps = 1000; // Real increment = 0.005; const UInt spatial_dimension = 3; Math::setTolerance(1.e-12); ElementType type = _tetrahedron_10; ElementType type_facet = Mesh::getFacetType(type); ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); Mesh mesh(spatial_dimension); mesh.read("tetrahedron.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull( SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); const MaterialCohesiveLinear<3> & mat_cohesive = dynamic_cast &>(model.getMaterial(1)); std::cout << mat_cohesive << std::endl; std::cout << model.getMaterial(2) << std::endl; const Real sigma_c = mat_cohesive.get("sigma_c"); const Real beta = mat_cohesive.get("beta"); + std::cout << sigma_c << " " << beta << std::endl; + Array & position = mesh.getNodes(); /* ------------------------------------------------------------------------ */ /* Facet part */ /* ------------------------------------------------------------------------ */ /// compute quadrature points positions on facets const Mesh & mesh_facets = model.getMeshFacets(); UInt nb_facet = mesh_facets.getNbElement(type_facet); UInt nb_quad_per_facet = model.getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); UInt nb_tot_quad = nb_quad_per_facet * nb_facet; Array quad_facets(nb_tot_quad, spatial_dimension); model.getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets, spatial_dimension, type_facet); /* ------------------------------------------------------------------------ */ /* End of facet part */ /* ------------------------------------------------------------------------ */ /// compute quadrature points position of the elements UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); UInt nb_element = mesh.getNbElement(type); UInt nb_tot_quad_el = nb_quad_per_element * nb_element; Array quad_elements(nb_tot_quad_el, spatial_dimension); model.getFEEngine().interpolateOnIntegrationPoints(position, quad_elements, spatial_dimension, type); /// assign some values to stresses Array & stress = const_cast &>(model.getMaterial(0).getStress(type)); Array::iterator> stress_it = stress.begin(spatial_dimension, spatial_dimension); for (UInt q = 0; q < nb_tot_quad_el; ++q, ++stress_it) { /// compute values for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = i; j < spatial_dimension; ++j) { UInt index = i * spatial_dimension + j; (*stress_it)(i, j) = index * function(sigma_c * 5, quad_elements(q, 0), quad_elements(q, 1), quad_elements(q, 2)); } } /// fill symmetrical part for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < i; ++j) { (*stress_it)(i, j) = (*stress_it)(j, i); } } } /// compute stress on facet quads Array stress_facets(nb_tot_quad, spatial_dimension * spatial_dimension); Array::iterator> stress_facets_it = stress_facets.begin(spatial_dimension, spatial_dimension); for (UInt q = 0; q < nb_tot_quad; ++q, ++stress_facets_it) { /// compute values for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = i; j < spatial_dimension; ++j) { UInt index = i * spatial_dimension + j; (*stress_facets_it)(i, j) = index * function(sigma_c * 5, quad_facets(q, 0), quad_facets(q, 1), quad_facets(q, 2)); } } /// fill symmetrical part for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < i; ++j) { (*stress_facets_it)(i, j) = (*stress_facets_it)(j, i); } } } /// insert cohesive elements model.checkCohesiveStress(); /// check insertion stress const Array & normals = model.getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(type_facet); const Array & tangents = model.getTangents(type_facet); const Array & sigma_c_eff = mat_cohesive.getInsertionTraction(type_cohesive); Vector normal_stress(spatial_dimension); const Array> & coh_element_to_facet = mesh_facets.getElementToSubelement(type_facet); Array::iterator> quad_facet_stress = stress_facets.begin(spatial_dimension, spatial_dimension); Array::const_iterator> quad_normal = normals.begin(spatial_dimension); Array::const_iterator> quad_tangents = tangents.begin(tangents.getNbComponent()); for (UInt f = 0; f < nb_facet; ++f) { const Element & cohesive_element = coh_element_to_facet(f)[1]; for (UInt q = 0; q < nb_quad_per_facet; ++q, ++quad_facet_stress, ++quad_normal, ++quad_tangents) { if (cohesive_element == ElementNull) continue; normal_stress.mul(*quad_facet_stress, *quad_normal); Real normal_contrib = normal_stress.dot(*quad_normal); Real first_tangent_contrib = 0; for (UInt dim = 0; dim < spatial_dimension; ++dim) first_tangent_contrib += normal_stress(dim) * (*quad_tangents)(dim); Real second_tangent_contrib = 0; for (UInt dim = 0; dim < spatial_dimension; ++dim) second_tangent_contrib += normal_stress(dim) * (*quad_tangents)(dim + spatial_dimension); Real tangent_contrib = std::sqrt(first_tangent_contrib * first_tangent_contrib + second_tangent_contrib * second_tangent_contrib); normal_contrib = std::max(0., normal_contrib); Real effective_norm = std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib / beta / beta); if (effective_norm < sigma_c) continue; if (!Math::are_float_equal( effective_norm, sigma_c_eff(cohesive_element.element * nb_quad_per_facet + q))) { std::cout << "Insertion tractions do not match" << std::endl; finalize(); return EXIT_FAILURE; } } } finalize(); std::cout << "OK: test_cohesive_extrinsic was passed!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat index bbd74ab72..4f6fc41e8 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/input_file.dat @@ -1,67 +1,65 @@ material elastic [ name = bulk rho = 2500 nu = 0.29 E = 70e9 # finite_deformation = 1 ] material cohesive_exponential [ name = coh1 sigma_c = 1.5e6 beta = 1 delta_c = 1e-4 exponential_penalty = true contact_tangent = 1.0 ] material cohesive_exponential [ name = coh2 sigma_c = 1.5e6 beta = 1 delta_c = 1e-4 exponential_penalty = true contact_tangent = 1.0 ] material cohesive_exponential [ name = coh3 sigma_c = 1.5e6 beta = 1 delta_c = 1e-4 exponential_penalty = true contact_tangent = 1.0 ] material cohesive_exponential [ name = coh4 sigma_c = 1.5e6 beta = 1 delta_c = 1e-4 exponential_penalty = true contact_tangent = 1.0 ] material cohesive_exponential [ name = coh5 sigma_c = 1.5e6 beta = 1 delta_c = 1e-4 exponential_penalty = true contact_tangent = 1.0 ] material cohesive_exponential [ name = interface sigma_c = 1.5e6 beta = 1 delta_c = 1e-4 exponential_penalty = true contact_tangent = 1.0 ] mesh parameters [ - cohesive_surfaces = coh1,coh2,coh3,coh4,coh5,interface - -] \ No newline at end of file +] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc index 02e268059..7dd425c39 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_insertion/test_cohesive_insertion_along_physical_surfaces.cc @@ -1,99 +1,107 @@ /** * @file test_cohesive_insertion_along_physical_surfaces.cc * * @author Fabian Barras * * @date creation: Fri Aug 07 2015 * * @brief Test intrinsic insertion of cohesive elements along physical surfaces * * @section LICENSE * * Copyright (©) 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 #include #include +#include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" +#include "material.hh" +#include "material_cohesive.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model_cohesive.hh" -#include "material.hh" -#include "material_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; -int main(int argc, char *argv[]) { +int main(int argc, char * argv[]) { initialize("input_file.dat", argc, argv); Math::setTolerance(1e-15); - + const UInt spatial_dimension = 3; Mesh mesh(spatial_dimension); mesh.read("3d_spherical_inclusion.msh"); - + SolidMechanicsModelCohesive model(mesh); - mesh.createGroupsFromMeshData("physical_names"); - model.initFull(SolidMechanicsModelCohesiveOptions(_static)); - - std::vector surfaces_name = {"interface", "coh1", "coh2", "coh3", "coh4", "coh5"}; + model.initFull(_analysis_method = _static); + + std::vector surfaces_name = {"interface", "coh1", "coh2", + "coh3", "coh4", "coh5"}; UInt nb_surf = surfaces_name.size(); - Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); - Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); + Mesh::type_iterator it = + mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); + Mesh::type_iterator end = + mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); - for(; it != end; ++it) { + for (; it != end; ++it) { for (UInt i = 0; i < nb_surf; ++i) { - - UInt expected_insertion = mesh.getElementGroup(surfaces_name[i]).getElements(mesh.getFacetType(*it)).size(); - UInt inserted_elements = model.getMaterial(surfaces_name[i]).getElementFilter()(*it).size(); + + UInt expected_insertion = mesh.getElementGroup(surfaces_name[i]) + .getElements(mesh.getFacetType(*it)) + .size(); + UInt inserted_elements = + model.getMaterial(surfaces_name[i]).getElementFilter()(*it).size(); AKANTU_DEBUG_ASSERT((expected_insertion == inserted_elements), - std::endl << "!!! Mismatch in insertion of surface named " - << surfaces_name[i] - << " --> " << inserted_elements << " inserted elements out of " - << expected_insertion << std::endl); - } + std::endl + << "!!! Mismatch in insertion of surface named " + << surfaces_name[i] << " --> " + << inserted_elements + << " inserted elements out of " + << expected_insertion << std::endl); + } } - /*std::string paraview_folder = "paraview/test_intrinsic_insertion_along_physical_surfaces/"; + /*std::string paraview_folder = + "paraview/test_intrinsic_insertion_along_physical_surfaces/"; model.setDirectory(paraview_folder); model.setBaseName("bulk"); model.addDumpField("partitions"); model.dump(); model.setDirectoryToDumper("cohesive elements", paraview_folder); model.setBaseNameToDumper("cohesive elements", "one_cohesive_element"); model.addDumpFieldToDumper("cohesive elements", "partitions"); model.addDumpFieldToDumper("cohesive elements", "material_index"); model.dump("cohesive elements"); */ model.assembleStiffnessMatrix(); finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc index b43a04742..dc4b498f2 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc @@ -1,430 +1,351 @@ /** * @file test_cohesive_intrinsic_tetrahedron.cc * * @author Marco Vocialta * * @date creation: Tue Aug 27 2013 * @date last modification: Thu Oct 15 2015 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; -void updateDisplacement(SolidMechanicsModelCohesive &, Array &, - ElementType, Vector &); +class Checker { +public: + Checker(const SolidMechanicsModelCohesive & model, + const Array & elements, ElementType type); + + void check(const Vector & opening, const Matrix & rotation) { + checkTractions(opening, rotation); + checkEquilibrium(); + computeEnergy(opening); + } + + void updateDisplacement(const Vector & increment); + +protected: + void checkTractions(const Vector & opening, const Matrix & rotation); + void checkEquilibrium(); + void checkResidual(const Matrix & rotation); + void computeEnergy(const Vector & opening); +private: + std::set nodes_to_check; + const SolidMechanicsModelCohesive & model; + ElementType type; + const Array & elements; + + const Material & mat_cohesive; + + Real sigma_c; + const Real beta; + const Real G_c; + const Real delta_0; + const Real kappa; + Real delta_c; + + const UInt spatial_dimension; -bool checkTractions(SolidMechanicsModelCohesive & model, ElementType type, - Vector & opening, Vector & theoretical_traction, - Matrix & rotation); + const ElementType type_facet; + const ElementType type_cohesive; + const Array & traction; + const Array & damage; + const UInt nb_quad_per_el; + const UInt nb_element; -void findNodesToCheck(const Mesh & mesh, const Array & elements, - ElementType type, Array & nodes_to_check); + const Real beta2_kappa2; + const Real beta2_kappa; -bool checkEquilibrium(const Array & residual); + Vector theoretical_traction; + Vector traction_old; + Vector opening_old; -bool checkResidual(const Array & residual, const Vector & traction, - const Array & nodes_to_check, - const Matrix & rotation); + Real Ed; +}; +/* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize("material_tetrahedron.dat", argc, argv); // debug::setDebugLevel(dblDump); - const UInt spatial_dimension = 3; const UInt max_steps = 60; const Real increment_constant = 0.01; Math::setTolerance(1.e-12); const ElementType type = _tetrahedron_10; Mesh mesh(spatial_dimension); mesh.read("tetrahedron.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(); - model.limitInsertion(_x, -0.01, 0.01); model.insertIntrinsicElements(); Array & boundary = model.getBlockedDOFs(); boundary.set(true); - UInt nb_element = mesh.getNbElement(type); - model.assembleInternalForces(); - model.setBaseName("intrinsic_tetrahedron"); model.addDumpFieldVector("displacement"); model.addDumpField("internal_force"); model.dump(); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_tetrahedron"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump("cohesive elements"); /// find elements to displace Array elements; - Real * bary = new Real[spatial_dimension]; + Vector bary(spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { - mesh.getBarycenter(el, type, bary); - if (bary[0] > 0.01) + mesh.getBarycenter(el, type, bary.storage()); + if (bary(_x) > 0.01) elements.push_back(el); } - delete[] bary; - - /// find nodes to check - Array nodes_to_check; - findNodesToCheck(mesh, elements, type, nodes_to_check); /// rotate mesh Real angle = 1.; - Matrix rotation(spatial_dimension, spatial_dimension); - rotation.clear(); - rotation(0, 0) = std::cos(angle); - rotation(0, 1) = std::sin(angle) * -1.; - rotation(1, 0) = std::sin(angle); - rotation(1, 1) = std::cos(angle); - rotation(2, 2) = 1.; - - Vector increment_tmp(spatial_dimension); - for (UInt dim = 0; dim < spatial_dimension; ++dim) { - increment_tmp(dim) = (dim + 1) * increment_constant; - } + // clang-format off + Matrix rotation{ + {std::cos(angle), std::sin(angle) * -1., 0.}, + {std::sin(angle), std::cos(angle), 0.}, + {0., 0., 1.}}; + // clang-format on - Vector increment(spatial_dimension); - increment.mul(rotation, increment_tmp); + Vector increment_tmp{increment_constant, 2. * increment_constant, + 3. * increment_constant}; + Vector increment = rotation * increment_tmp; - Array & position = mesh.getNodes(); - Array position_tmp(position); + auto & position = mesh.getNodes(); + auto position_it = position.begin(spatial_dimension); + auto position_end = position.end(spatial_dimension); - Array::iterator> position_it = - position.begin(spatial_dimension); - Array::iterator> position_end = - position.end(spatial_dimension); - Array::iterator> position_tmp_it = - position_tmp.begin(spatial_dimension); - - for (; position_it != position_end; ++position_it, ++position_tmp_it) - position_it->mul(rotation, *position_tmp_it); + for (; position_it != position_end; ++position_it) { + auto & pos = *position_it; + pos = rotation * pos; + } model.dump(); model.dump("cohesive elements"); - updateDisplacement(model, elements, type, increment); - - Real theoretical_Ed = 0; - - Vector opening(spatial_dimension); - Vector traction(spatial_dimension); - Vector opening_old(spatial_dimension); - Vector traction_old(spatial_dimension); + /// find nodes to check + Checker checker(model, elements, type); - opening.clear(); - traction.clear(); - opening_old.clear(); - traction_old.clear(); + checker.updateDisplacement(increment); - Vector Dt(spatial_dimension); - Vector Do(spatial_dimension); + Real theoretical_Ed = 0; - const Array & residual = model.getInternalForce(); + Vector opening(spatial_dimension, 0.); + Vector opening_old(spatial_dimension, 0.); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { + model.solveStep(); - model.assembleInternalForces(); + model.dump(); + model.dump("cohesive elements"); opening += increment_tmp; - if (checkTractions(model, type, opening, traction, rotation) || - checkEquilibrium(residual) || - checkResidual(residual, traction, nodes_to_check, rotation)) { - finalize(); - return EXIT_FAILURE; - } - - /// compute energy - Do = opening; - Do -= opening_old; - - Dt = traction_old; - Dt += traction; - - theoretical_Ed += .5 * Do.dot(Dt); - - opening_old = opening; - traction_old = traction; + checker.check(opening, rotation); - updateDisplacement(model, elements, type, increment); - - if (s % 10 == 0) { - std::cout << "passing step " << s << "/" << max_steps << std::endl; - model.dump(); - model.dump("cohesive elements"); - } + checker.updateDisplacement(increment); } model.dump(); model.dump("cohesive elements"); Real Ed = model.getEnergy("dissipated"); theoretical_Ed *= 4.; std::cout << Ed << " " << theoretical_Ed << std::endl; if (!Math::are_float_equal(Ed, theoretical_Ed) || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic_tetrahedron was passed!" << std::endl; return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ - -void updateDisplacement(SolidMechanicsModelCohesive & model, - Array & elements, ElementType type, - Vector & increment) { - - UInt spatial_dimension = model.getSpatialDimension(); +void Checker::updateDisplacement(const Vector & increment) { Mesh & mesh = model.getFEEngine().getMesh(); - UInt nb_element = elements.size(); - UInt nb_nodes = mesh.getNbNodes(); - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - - const Array & connectivity = mesh.getConnectivity(type); - Array & displacement = model.getDisplacement(); - Array update(nb_nodes); + const auto & connectivity = mesh.getConnectivity(type); + auto & displacement = model.getDisplacement(); + Array update(displacement.size()); update.clear(); - for (UInt el = 0; el < nb_element; ++el) { - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt node = connectivity(elements(el), n); + auto conn_it = connectivity.begin(connectivity.getNbComponent()); + auto conn_end = connectivity.begin(connectivity.getNbComponent()); + for (;conn_it != conn_end;++conn_it) { + const auto & conn = *conn_it; + for (UInt n = 0; n < conn.size(); ++n) { + UInt node = conn(n); if (!update(node)) { Vector node_disp(displacement.storage() + node * spatial_dimension, spatial_dimension); node_disp += increment; update(node) = true; } } } } /* -------------------------------------------------------------------------- */ +Checker::Checker(const SolidMechanicsModelCohesive & model, + const Array & elements, ElementType type) + : model(model), type(std::move(type)), elements(elements), + mat_cohesive(model.getMaterial(1)), sigma_c(mat_cohesive.get("sigma_c")), + beta(mat_cohesive.get("beta")), G_c(mat_cohesive.get("G_c")), + delta_0(mat_cohesive.get("delta_0")), kappa(mat_cohesive.get("kappa")), + spatial_dimension(model.getSpatialDimension()), + type_facet(Mesh::getFacetType(type)), + type_cohesive(FEEngine::getCohesiveElementType(type_facet)), + traction(mat_cohesive.getArray("tractions", type_cohesive)), + damage(mat_cohesive.getArray("damage", type_cohesive)), + nb_quad_per_el(model.getFEEngine("CohesiveFEEngine") + .getNbIntegrationPoints(type_cohesive)), + nb_element(model.getMesh().getNbElement(type_cohesive)), + beta2_kappa2(beta * beta / kappa / kappa), + beta2_kappa(beta * beta / kappa) { + const Mesh & mesh = model.getMesh(); + const auto & connectivity = mesh.getConnectivity(type); + const auto & position = mesh.getNodes(); + auto conn_it = connectivity.begin(connectivity.getNbComponent()); + for (const auto & element : elements) { + Vector conn_el(conn_it[element]); + for (UInt n = 0; n < conn_el.size(); ++n) { + UInt node = conn_el(n); + if (Math::are_float_equal(position(node, _x), 0.)) + nodes_to_check.insert(node); + } + } -bool checkTractions(SolidMechanicsModelCohesive & model, ElementType type, - Vector & opening, Vector & theoretical_traction, - Matrix & rotation) { - UInt spatial_dimension = model.getSpatialDimension(); - - const MaterialCohesive & mat_cohesive = - dynamic_cast(model.getMaterial(1)); - - Real sigma_c = mat_cohesive.getParam("sigma_c"); - const Real beta = mat_cohesive.getParam("beta"); - const Real G_c = mat_cohesive.getParam("G_c"); - const Real delta_0 = mat_cohesive.getParam("delta_0"); - const Real kappa = mat_cohesive.getParam("kappa"); - Real delta_c = 2 * G_c / sigma_c; + delta_c = 2 * G_c / sigma_c; sigma_c *= delta_c / (delta_c - delta_0); +} - ElementType type_facet = Mesh::getFacetType(type); - ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); - - const Array & traction = mat_cohesive.getTraction(type_cohesive); - const Array & damage = mat_cohesive.getDamage(type_cohesive); - - UInt nb_quad_per_el = model.getFEEngine("CohesiveFEEngine") - .getNbIntegrationPoints(type_cohesive); - UInt nb_element = model.getMesh().getNbElement(type_cohesive); - UInt tot_nb_quad = nb_element * nb_quad_per_el; - - Vector normal_opening(spatial_dimension); - normal_opening.clear(); - normal_opening(0) = opening(0); - Real normal_opening_norm = normal_opening.norm(); - - Vector tangential_opening(spatial_dimension); - tangential_opening.clear(); - for (UInt dim = 1; dim < spatial_dimension; ++dim) - tangential_opening(dim) = opening(dim); - - Real tangential_opening_norm = tangential_opening.norm(); - - Real beta2_kappa2 = beta * beta / kappa / kappa; - Real beta2_kappa = beta * beta / kappa; - - Real delta = std::sqrt(tangential_opening_norm * tangential_opening_norm * +/* -------------------------------------------------------------------------- */ +void Checker::checkTractions(const Vector & opening, const Matrix & rotation) { + auto normal_opening = opening * Vector{1., 0., 0.}; + auto tangential_opening = opening - normal_opening; + const Real normal_opening_norm = normal_opening.norm(); + const Real tangential_opening_norm = tangential_opening.norm(); + const Real delta = + std::max(std::sqrt(tangential_opening_norm * tangential_opening_norm * beta2_kappa2 + - normal_opening_norm * normal_opening_norm); - - delta = std::max(delta, delta_0); + normal_opening_norm * normal_opening_norm), + 0.); Real theoretical_damage = std::min(delta / delta_c, 1.); - if (Math::are_float_equal(theoretical_damage, 1.)) - theoretical_traction.clear(); - else { - theoretical_traction = tangential_opening; - theoretical_traction *= beta2_kappa; - theoretical_traction += normal_opening; - theoretical_traction *= sigma_c / delta * (1. - theoretical_damage); - } - + theoretical_traction = (tangential_opening * beta2_kappa + normal_opening) * + sigma_c / delta * (1. - theoretical_damage); // adjust damage theoretical_damage = std::max((delta - delta_0) / (delta_c - delta_0), 0.); theoretical_damage = std::min(theoretical_damage, 1.); - Vector theoretical_traction_rotated(spatial_dimension); - theoretical_traction_rotated.mul(rotation, theoretical_traction); - - for (UInt q = 0; q < tot_nb_quad; ++q) { - for (UInt dim = 0; dim < spatial_dimension; ++dim) { - if (!Math::are_float_equal(theoretical_traction_rotated(dim), - traction(q, dim))) { - std::cout << "Tractions are incorrect" << std::endl; - return 1; - } - } - - if (!Math::are_float_equal(theoretical_damage, damage(q))) { - std::cout << "Damage is incorrect" << std::endl; - return 1; - } - } - - return 0; + Vector theoretical_traction_rotated = rotation * theoretical_traction; + + std::for_each( + traction.begin(spatial_dimension), traction.end(spatial_dimension), + [&theoretical_traction_rotated](auto && traction) { + Real diff = Vector(theoretical_traction_rotated - traction).norm(); + if (diff > 1e-14) + throw std::domain_error("Tractions are incorrect"); + }); + + std::for_each(damage.begin(), damage.end(), + [&theoretical_damage](auto && damage) { + if (not Math::are_float_equal(theoretical_damage, damage)) + throw std::domain_error("Damage is incorrect"); + }); } /* -------------------------------------------------------------------------- */ +void Checker::computeEnergy(const Vector & opening) { + /// compute energy + auto Do = opening - opening_old; + auto Dt = traction_old + theoretical_traction; -void findNodesToCheck(const Mesh & mesh, const Array & elements, - ElementType type, Array & nodes_to_check) { - - const Array & connectivity = mesh.getConnectivity(type); - const Array & position = mesh.getNodes(); - - UInt nb_nodes = position.size(); - UInt nb_nodes_per_elem = connectivity.getNbComponent(); - - Array checked_nodes(nb_nodes); - checked_nodes.clear(); - - for (UInt el = 0; el < elements.size(); ++el) { - - UInt element = elements(el); - Vector conn_el(connectivity.storage() + nb_nodes_per_elem * element, - nb_nodes_per_elem); + Ed += .5 * Do.dot(Dt); - for (UInt n = 0; n < nb_nodes_per_elem; ++n) { - UInt node = conn_el(n); - if (Math::are_float_equal(position(node, 0), 0.) && - checked_nodes(node) == false) { - checked_nodes(node) = true; - nodes_to_check.push_back(node); - } - } - } + opening_old = opening; + traction_old = theoretical_traction; } /* -------------------------------------------------------------------------- */ - -bool checkEquilibrium(const Array & residual) { - - UInt spatial_dimension = residual.getNbComponent(); - - Vector residual_sum(spatial_dimension); - residual_sum.clear(); - - Array::const_iterator> res_it = - residual.begin(spatial_dimension); - Array::const_iterator> res_end = - residual.end(spatial_dimension); - +void Checker::checkEquilibrium() { + Vector residual_sum(spatial_dimension, 0.); + const auto & residual = model.getInternalForce(); + auto res_it = residual.begin(spatial_dimension); + auto res_end = residual.end(spatial_dimension); for (; res_it != res_end; ++res_it) residual_sum += *res_it; - for (UInt s = 0; s < spatial_dimension; ++s) { - if (!Math::are_float_equal(residual_sum(s), 0.)) { - std::cout << "System is not in equilibrium!" << std::endl; - return 1; - } - } - - return 0; + if (!Math::are_float_equal(residual_sum.norm(), 0.)) + throw std::domain_error("System is not in equilibrium!"); } /* -------------------------------------------------------------------------- */ - -bool checkResidual(const Array & residual, const Vector & traction, - const Array & nodes_to_check, - const Matrix & rotation) { - - UInt spatial_dimension = residual.getNbComponent(); - - Vector total_force(spatial_dimension); - total_force.clear(); - - for (UInt n = 0; n < nodes_to_check.size(); ++n) { - UInt node = nodes_to_check(n); - - Vector res(residual.storage() + node * spatial_dimension, - spatial_dimension); - +void Checker::checkResidual(const Matrix & rotation) { + Vector total_force(spatial_dimension, 0.); + const auto & residual = model.getInternalForce(); + for (auto node : nodes_to_check) { + Vector res(residual.begin(spatial_dimension)[node]); total_force += res; } Vector theoretical_total_force(spatial_dimension); - theoretical_total_force.mul(rotation, traction); + theoretical_total_force.mul(rotation, theoretical_traction); theoretical_total_force *= -1 * 2 * 2; for (UInt s = 0; s < spatial_dimension; ++s) { if (!Math::are_float_equal(total_force(s), theoretical_total_force(s))) { std::cout << "Total force isn't correct!" << std::endl; - return 1; + std::terminate(); } } - - return 0; }