diff --git a/src/fe_engine/element_class.hh b/src/fe_engine/element_class.hh index e8eac82ac..fce855a3f 100644 --- a/src/fe_engine/element_class.hh +++ b/src/fe_engine/element_class.hh @@ -1,398 +1,399 @@ /** * @file element_class.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 21 2016 * * @brief Declaration of the ElementClass main class and the * Integration and Interpolation elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_HH__ #define __AKANTU_ELEMENT_CLASS_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /// default element class structure template struct ElementClassProperty { static const GeometricalType geometrical_type = _gt_not_defined; static const InterpolationType interpolation_type = _itp_not_defined; static const ElementKind element_kind = _ek_regular; static const UInt spatial_dimension = 0; static const GaussIntegrationType gauss_integration_type = _git_not_defined; static const UInt polynomial_degree = 0; }; /// Macro to generate the element class structures for different element types #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type = geom_type; \ static const InterpolationType interpolation_type = interp_type; \ static const ElementKind element_kind = elem_kind; \ static const UInt spatial_dimension = sp; \ static const GaussIntegrationType gauss_integration_type = \ gauss_int_type; \ static const UInt polynomial_degree = min_int_order; \ } /* -------------------------------------------------------------------------- */ /* Geometry */ /* -------------------------------------------------------------------------- */ /// Default GeometricalShape structure template struct GeometricalShape { static const GeometricalShapeType shape = _gst_point; }; /// Templated GeometricalShape with function contains template struct GeometricalShapeContains { /// Check if the point (vector in 2 and 3D) at natural coordinate coor template static inline bool contains(const vector_type & coord); }; /// Macro to generate the GeometricalShape structures for different geometrical /// types #define AKANTU_DEFINE_SHAPE(geom_type, geom_shape) \ template <> struct GeometricalShape { \ static const GeometricalShapeType shape = geom_shape; \ } /* -------------------------------------------------------------------------- */ /// Templated GeometricalElement with function getInradius template ::shape> class GeometricalElement { public: /// compute the in-radius static inline Real getInradius(__attribute__((unused)) const Matrix & coord) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// true if the natural coordinates are in the element template static inline bool contains(const vector_type & coord); public: static AKANTU_GET_MACRO_NOT_CONST(SpatialDimension, spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbNodesPerElement, nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbFacetTypes, nb_facet_types, UInt); static inline UInt getNbFacetsPerElement(UInt t); static inline UInt getNbFacetsPerElement(); static inline const MatrixProxy getFacetLocalConnectivityPerElement(UInt t = 0); protected: /// Number of nodes per element static UInt nb_nodes_per_element; /// spatial dimension of the element static UInt spatial_dimension; /// number of different facet types static UInt nb_facet_types; /// number of facets for element static UInt nb_facets[]; /// storage of the facet local connectivity static UInt facet_connectivity_vect[]; /// local connectivity of facets static UInt * facet_connectivity[]; private: /// Type of the facet elements static UInt nb_nodes_per_facet[]; }; /* -------------------------------------------------------------------------- */ /* Interpolation */ /* -------------------------------------------------------------------------- */ /// default InterpolationPorperty structure template struct InterpolationPorperty { static const InterpolationKind kind = _itk_not_defined; static const UInt nb_nodes_per_element = 0; static const UInt natural_space_dimension = 0; }; /// Macro to generate the InterpolationPorperty structures for different /// interpolation types #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) \ template <> struct InterpolationPorperty { \ static const InterpolationKind kind = itp_kind; \ static const UInt nb_nodes_per_element = nb_nodes; \ static const UInt natural_space_dimension = ndim; \ } #include "interpolation_element_tmpl.hh" /* -------------------------------------------------------------------------- */ /// Generic (templated by the enum InterpolationType which specifies the order /// and the dimension of the interpolation) class handling the elemental /// interpolation template ::kind> class InterpolationElement { public: typedef InterpolationPorperty interpolation_property; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix & natural_coord, Matrix & N); /// compute the shape values for a given point in natural coordinates template static inline void computeShapes(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) vector_type & N) { AKANTU_DEBUG_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix & natural_coord, Tensor3 & dnds); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template static inline void computeDNDS(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) matrix_type & dnds) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute jacobian (or integration variable change factor) for a given point /// in the case of spatial_dimension != natural_space_dimension static inline void computeSpecialJacobian(__attribute__((unused)) const Matrix & J, __attribute__((unused)) Real & jacobians) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate a field given (arbitrary) natural coordinates static inline void interpolateOnNaturalCoordinates(const Vector & natural_coords, const Matrix & nodal_values, Vector & interpolated); /// interpolate a field given the shape functions on the interpolation point static inline void interpolate(const Matrix & nodal_values, const Vector & shapes, Vector & interpolated); /// interpolate a field given the shape functions on the interpolations points static inline void interpolate(const Matrix & nodal_values, const Matrix & shapes, Matrix & interpolated); /// compute the gradient of a given field on the given natural coordinates static inline void gradientOnNaturalCoordinates(const Vector & natural_coords, const Matrix & f, Matrix & gradient); public: static AKANTU_GET_MACRO_NOT_CONST( ShapeSize, InterpolationPorperty::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeDerivativesSize, (InterpolationPorperty::nb_nodes_per_element * InterpolationPorperty::natural_space_dimension), UInt); static AKANTU_GET_MACRO_NOT_CONST( NaturalSpaceDimension, InterpolationPorperty::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbNodesPerInterpolationElement, InterpolationPorperty::nb_nodes_per_element, UInt); }; /* -------------------------------------------------------------------------- */ /* Integration */ /* -------------------------------------------------------------------------- */ template struct GaussIntegrationTypeData { /// quadrature points in natural coordinates static Real quad_positions[]; /// weights for the Gauss integration static Real quad_weights[]; }; template ::polynomial_degree> class GaussIntegrationElement { public: static UInt getNbQuadraturePoints(); static const Matrix getQuadraturePoints(); static const Vector getWeights(); }; /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ template ::element_kind> class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: typedef GeometricalElement< ElementClassProperty::geometrical_type> geometrical_element; typedef InterpolationElement::interpolation_type> interpolation_element; typedef ElementClassProperty element_property; typedef typename interpolation_element::interpolation_property interpolation_property; public: /** * compute @f$ J = \frac{\partial x_j}{\partial s_i} @f$ the variation of real * coordinates along with variation of natural coordinates on a given point in * natural coordinates */ static inline void computeJMat(const Matrix & dnds, const Matrix & node_coords, Matrix & J); /** * compute the Jacobian matrix by computing the variation of real coordinates * along with variation of natural coordinates on a given set of points in * natural coordinates */ static inline void computeJMat(const Tensor3 & dnds, const Matrix & node_coords, Tensor3 & J); /// compute the jacobians of a serie of natural coordinates static inline void computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians); /// compute jacobian (or integration variable change factor) for a set of /// points static inline void computeJacobian(const Tensor3 & J, Vector & jacobians); /// compute jacobian (or integration variable change factor) for a given point static inline void computeJacobian(const Matrix & J, Real & jacobians); /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3 & J, const Tensor3 & dnds, Tensor3 & shape_deriv); /// compute shape derivatives (input is dxds) for a given point static inline void computeShapeDerivatives(const Matrix & J, const Matrix & dnds, Matrix & shape_deriv); /// compute the normal of a surface defined by the function f static inline void computeNormalsOnNaturalCoordinates(const Matrix & coord, Matrix & f, Matrix & normals); /// get natural coordinates from real coordinates static inline void inverseMap(const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, Real tolerance = 1e-8); /// get natural coordinates from real coordinates static inline void inverseMap(const Matrix & real_coords, const Matrix & node_coords, Matrix & natural_coords, Real tolerance = 1e-8); public: static AKANTU_GET_MACRO_NOT_CONST(Kind, element_kind, ElementKind); static AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty::spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, p1_type, const ElementType &); static const ElementType & getFacetType(UInt t = 0) { return facet_type[t]; } static ElementType * getFacetTypeInternal() { return facet_type; } protected: /// Type of the facet elements static ElementType facet_type[]; /// type of element P1 associated static ElementType p1_type; }; /* -------------------------------------------------------------------------- */ } // akantu #include "element_class_tmpl.hh" - - /* -------------------------------------------------------------------------- */ namespace akantu { #include "element_class_point_1_inline_impl.cc" #include "element_class_segment_2_inline_impl.cc" #include "element_class_segment_3_inline_impl.cc" #include "element_class_triangle_3_inline_impl.cc" #include "element_class_triangle_6_inline_impl.cc" #include "element_class_tetrahedron_4_inline_impl.cc" #include "element_class_tetrahedron_10_inline_impl.cc" #include "element_class_quadrangle_4_inline_impl.cc" #include "element_class_quadrangle_8_inline_impl.cc" #include "element_class_hexahedron_8_inline_impl.cc" #include "element_class_hexahedron_20_inline_impl.cc" #include "element_class_pentahedron_6_inline_impl.cc" #include "element_class_pentahedron_15_inline_impl.cc" } // akantu /* -------------------------------------------------------------------------- */ +#if defined(AKANTU_COHESIVE_ELEMENT) +#include "cohesive_element.hh" +#endif #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "element_class_structural.hh" #endif #if defined(AKANTU_IGFEM) #include "element_class_igfem.hh" #endif #endif /* __AKANTU_ELEMENT_CLASS_HH__ */ diff --git a/src/model/boundary_condition_tmpl.hh b/src/model/boundary_condition_tmpl.hh index 4957dd61e..8d42ed878 100644 --- a/src/model/boundary_condition_tmpl.hh +++ b/src/model/boundary_condition_tmpl.hh @@ -1,270 +1,268 @@ /** * @file boundary_condition_tmpl.hh * * @author Dana Christen * @author Nicolas Richart * * @date creation: Fri May 03 2013 * @date last modification: Fri Oct 16 2015 * * @brief implementation of the applyBC * * @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 "element_group.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template void BoundaryCondition::initBC(ModelType & model, Array & primal, Array & dual) { this->model = &model; this->primal = &primal; this->dual = &dual; if (this->model->getSpatialDimension() > 1) this->model->initFEEngineBoundary(); } /* -------------------------------------------------------------------------- */ template void BoundaryCondition::initBC(ModelType & model, Array & primal, Array & primal_increment, Array & dual) { this->initBC(model, primal, dual); this->primal_increment = &primal_increment; } /* -------------------------------------------------------------------------- */ /* Partial specialization for DIRICHLET functors */ template template struct BoundaryCondition::TemplateFunctionWrapper< FunctorType, BC::Functor::_dirichlet> { static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance) { ModelType & model = bc_instance.getModel(); Array & primal = bc_instance.getPrimal(); const Array & coords = model.getMesh().getNodes(); Array & boundary_flags = model.getBlockedDOFs(); UInt dim = model.getMesh().getSpatialDimension(); Array::vector_iterator primal_iter = primal.begin(primal.getNbComponent()); Array::const_vector_iterator coords_iter = coords.begin(dim); Array::vector_iterator flags_iter = boundary_flags.begin(boundary_flags.getNbComponent()); - for (ElementGroup::const_node_iterator nodes_it(group.node_begin()); - nodes_it != group.node_end(); ++nodes_it) { - UInt n = *nodes_it; + for (auto n : group.getNodeGroup()) { Vector flag(flags_iter[n]); Vector primal(primal_iter[n]); Vector coords(coords_iter[n]); func(n, flag, primal, coords); } } }; /* -------------------------------------------------------------------------- */ /* Partial specialization for NEUMANN functors */ template template struct BoundaryCondition::TemplateFunctionWrapper< FunctorType, BC::Functor::_neumann> { static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance) { UInt dim = bc_instance.getModel().getSpatialDimension(); switch (dim) { case 1: { AKANTU_DEBUG_TO_IMPLEMENT(); break; } case 2: case 3: { applyBC(func, group, bc_instance, _not_ghost); applyBC(func, group, bc_instance, _ghost); break; } } } static inline void applyBC(const FunctorType & func, const ElementGroup & group, BoundaryCondition & bc_instance, GhostType ghost_type) { ModelType & model = bc_instance.getModel(); Array & dual = bc_instance.getDual(); const Mesh & mesh = model.getMesh(); const Array & nodes_coords = mesh.getNodes(); const FEEngine & fem_boundary = model.getFEEngineBoundary(); UInt dim = model.getSpatialDimension(); UInt nb_degree_of_freedom = dual.getNbComponent(); IntegrationPoint quad_point; quad_point.ghost_type = ghost_type; ElementGroup::type_iterator type_it = group.firstType(dim - 1, ghost_type); ElementGroup::type_iterator type_end = group.lastType(dim - 1, ghost_type); // Loop over the boundary element types for (; type_it != type_end; ++type_it) { const Array & element_ids = group.getElements(*type_it, ghost_type); Array::const_scalar_iterator elem_iter = element_ids.begin(); Array::const_scalar_iterator elem_iter_end = element_ids.end(); UInt nb_quad_points = fem_boundary.getNbIntegrationPoints(*type_it, ghost_type); UInt nb_elements = element_ids.getSize(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*type_it); Array * dual_before_integ = new Array( nb_elements * nb_quad_points, nb_degree_of_freedom, 0.); Array * quad_coords = new Array(nb_elements * nb_quad_points, dim); const Array & normals_on_quad = fem_boundary.getNormalsOnIntegrationPoints(*type_it, ghost_type); fem_boundary.interpolateOnIntegrationPoints( nodes_coords, *quad_coords, dim, *type_it, ghost_type, element_ids); Array::const_vector_iterator normals_begin = normals_on_quad.begin(dim); Array::const_vector_iterator normals_iter; Array::const_vector_iterator quad_coords_iter = quad_coords->begin(dim); Array::vector_iterator dual_iter = dual_before_integ->begin(nb_degree_of_freedom); quad_point.type = *type_it; for (; elem_iter != elem_iter_end; ++elem_iter) { UInt el = *elem_iter; quad_point.element = el; normals_iter = normals_begin + el * nb_quad_points; for (UInt q(0); q < nb_quad_points; ++q) { quad_point.num_point = q; func(quad_point, *dual_iter, *quad_coords_iter, *normals_iter); ++dual_iter; ++quad_coords_iter; ++normals_iter; } } delete quad_coords; /* -------------------------------------------------------------------- */ // Initialization of iterators Array::matrix_iterator dual_iter_mat = dual_before_integ->begin(nb_degree_of_freedom, 1); elem_iter = element_ids.begin(); Array::const_matrix_iterator shapes_iter_begin = fem_boundary.getShapes(*type_it, ghost_type) .begin(1, nb_nodes_per_element); Array * dual_by_shapes = new Array(nb_elements * nb_quad_points, nb_degree_of_freedom * nb_nodes_per_element); Array::matrix_iterator dual_by_shapes_iter = dual_by_shapes->begin(nb_degree_of_freedom, nb_nodes_per_element); Array::const_matrix_iterator shapes_iter; /* -------------------------------------------------------------------- */ // Loop computing dual x shapes for (; elem_iter != elem_iter_end; ++elem_iter) { shapes_iter = shapes_iter_begin + *elem_iter * nb_quad_points; for (UInt q(0); q < nb_quad_points; ++q, ++dual_iter_mat, ++dual_by_shapes_iter, ++shapes_iter) { dual_by_shapes_iter->mul(*dual_iter_mat, *shapes_iter); } } delete dual_before_integ; Array * dual_by_shapes_integ = new Array( nb_elements, nb_degree_of_freedom * nb_nodes_per_element); fem_boundary.integrate(*dual_by_shapes, *dual_by_shapes_integ, nb_degree_of_freedom * nb_nodes_per_element, *type_it, ghost_type, element_ids); delete dual_by_shapes; // assemble the result into force vector model.getDOFManager().assembleElementalArrayLocalArray( *dual_by_shapes_integ, dual, *type_it, ghost_type, 1., element_ids); delete dual_by_shapes_integ; } } }; /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func) { GroupManager::const_element_group_iterator bit = model->getMesh().getGroupManager().element_group_begin(); GroupManager::const_element_group_iterator bend = model->getMesh().getGroupManager().element_group_end(); for (; bit != bend; ++bit) applyBC(func, *bit); } /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func, const std::string & group_name) { try { const ElementGroup & element_group = model->getMesh().getElementGroup(group_name); applyBC(func, element_group); } catch (akantu::debug::Exception e) { AKANTU_EXCEPTION("Error applying a boundary condition onto \"" << group_name << "\"! [" << e.what() << "]"); } } /* -------------------------------------------------------------------------- */ template template inline void BoundaryCondition::applyBC(const FunctorType & func, const ElementGroup & element_group) { #if !defined(AKANTU_NDEBUG) if (element_group.getDimension() != model->getSpatialDimension() - 1) AKANTU_DEBUG_WARNING("The group " << element_group.getName() << " does not contain only boundaries elements"); #endif TemplateFunctionWrapper::applyBC(func, element_group, *this); } } // akantu diff --git a/src/model/common/neighborhood_base.cc b/src/model/common/neighborhood_base.cc index f43045e18..c695337c1 100644 --- a/src/model/common/neighborhood_base.cc +++ b/src/model/common/neighborhood_base.cc @@ -1,306 +1,305 @@ /** * @file neighborhood_base.cc * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Sat Sep 26 2015 * @date last modification: Wed Nov 25 2015 * * @brief Implementation of generic neighborhood base * * @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 "neighborhood_base.hh" #include "grid_synchronizer.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NeighborhoodBase::NeighborhoodBase(Model & model, const ElementTypeMapReal & quad_coordinates, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), model(model), neighborhood_radius(0.), spatial_grid(nullptr), is_creating_grid(false), grid_synchronizer(nullptr), quad_coordinates(quad_coordinates), spatial_dimension(this->model.getMesh().getSpatialDimension()) { AKANTU_DEBUG_IN(); this->registerDataAccessor(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NeighborhoodBase::~NeighborhoodBase() = default; /* -------------------------------------------------------------------------- */ // void NeighborhoodBase::createSynchronizerRegistry( // DataAccessor * data_accessor) { // this->synch_registry = new SynchronizerRegistry(*data_accessor); // } /* -------------------------------------------------------------------------- */ void NeighborhoodBase::initNeighborhood() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Creating the grid"); this->createGrid(); AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void NeighborhoodBase::createGrid() { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model.getMesh(); - mesh.computeBoundingBox(); - const Vector & lower_bounds = mesh.getLocalLowerBounds(); - const Vector & upper_bounds = mesh.getLocalUpperBounds(); + const auto & lower_bounds = mesh.getLocalLowerBounds(); + const auto & upper_bounds = mesh.getLocalUpperBounds(); Vector center = 0.5 * (upper_bounds + lower_bounds); Vector spacing(spatial_dimension, this->neighborhood_radius * safety_factor); spatial_grid = std::make_unique>(spatial_dimension, spacing, center); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NeighborhoodBase::updatePairList() { AKANTU_DEBUG_IN(); //// loop over all quads -> all cells SpatialGrid::cells_iterator cell_it = spatial_grid->beginCells(); SpatialGrid::cells_iterator cell_end = spatial_grid->endCells(); Vector q1_coords(spatial_dimension); Vector q2_coords(spatial_dimension); IntegrationPoint q1; IntegrationPoint q2; UInt counter = 0; for (; cell_it != cell_end; ++cell_it) { AKANTU_DEBUG_INFO("Looping on next cell"); SpatialGrid::Cell::iterator first_quad = spatial_grid->beginCell(*cell_it); SpatialGrid::Cell::iterator last_quad = spatial_grid->endCell(*cell_it); for (; first_quad != last_quad; ++first_quad, ++counter) { q1 = *first_quad; if (q1.ghost_type == _ghost) break; Array::const_vector_iterator coords_type_1_it = this->quad_coordinates(q1.type, q1.ghost_type) .begin(spatial_dimension); q1_coords = coords_type_1_it[q1.global_num]; AKANTU_DEBUG_INFO("Current quadrature point in this cell: " << q1); SpatialGrid::CellID cell_id = spatial_grid->getCellID(q1_coords); /// loop over all the neighbouring cells of the current quad SpatialGrid::neighbor_cells_iterator first_neigh_cell = spatial_grid->beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = spatial_grid->endNeighborCells(cell_id); for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { SpatialGrid::Cell::iterator first_neigh_quad = spatial_grid->beginCell(*first_neigh_cell); SpatialGrid::Cell::iterator last_neigh_quad = spatial_grid->endCell(*first_neigh_cell); // loop over the quadrature point in the current neighboring cell for (; first_neigh_quad != last_neigh_quad; ++first_neigh_quad) { q2 = *first_neigh_quad; Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type) .begin(spatial_dimension); q2_coords = coords_type_2_it[q2.global_num]; Real distance = q1_coords.distance(q2_coords); if (distance <= this->neighborhood_radius + Math::getTolerance() && (q2.ghost_type == _ghost || (q2.ghost_type == _not_ghost && q1.global_num <= q2.global_num))) { // storing only half lists pair_list[q2.ghost_type].push_back(std::make_pair(q1, q2)); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NeighborhoodBase::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType)gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); for (; first_pair != last_pair; ++first_pair) { const IntegrationPoint & q1 = first_pair->first; const IntegrationPoint & q2 = first_pair->second; pout << q1 << " " << q2 << " " << std::endl; } } } /* -------------------------------------------------------------------------- */ void NeighborhoodBase::saveNeighborCoords(const std::string & filename) const { /// this function is not optimazed and only used for tests on small meshes /// @todo maybe optimize this function for better performance? Vector q1_coords(spatial_dimension); Vector q2_coords(spatial_dimension); IntegrationPoint q1; IntegrationPoint q2; std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); /// loop over all the quads and write the position of their neighbors SpatialGrid::cells_iterator cell_it = spatial_grid->beginCells(); SpatialGrid::cells_iterator cell_end = spatial_grid->endCells(); for (; cell_it != cell_end; ++cell_it) { SpatialGrid::Cell::iterator first_quad = spatial_grid->beginCell(*cell_it); SpatialGrid::Cell::iterator last_quad = spatial_grid->endCell(*cell_it); for (; first_quad != last_quad; ++first_quad) { q1 = *first_quad; Array::const_vector_iterator coords_type_1_it = this->quad_coordinates(q1.type, q1.ghost_type) .begin(spatial_dimension); q1_coords = coords_type_1_it[q1.global_num]; pout << "#neighbors for quad " << q1.global_num << std::endl; pout << q1_coords << std::endl; for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType)gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); for (; first_pair != last_pair; ++first_pair) { if (q1 == first_pair->first && first_pair->second != q1) { q2 = first_pair->second; Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type) .begin(spatial_dimension); q2_coords = coords_type_2_it[q2.global_num]; pout << q2_coords << std::endl; } if (q1 == first_pair->second && first_pair->first != q1) { q2 = first_pair->first; Array::const_vector_iterator coords_type_2_it = this->quad_coordinates(q2.type, q2.ghost_type) .begin(spatial_dimension); q2_coords = coords_type_2_it[q2.global_num]; pout << q2_coords << std::endl; } } } } } } /* -------------------------------------------------------------------------- */ void NeighborhoodBase::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { AKANTU_DEBUG_IN(); FEEngine & fem = this->model.getFEEngine(); UInt nb_quad = 0; // Change the pairs in new global numbering for(auto ghost_type2 : ghost_types) { auto first_pair = pair_list[ghost_type2].begin(); auto last_pair = pair_list[ghost_type2].end(); for (; first_pair != last_pair; ++first_pair) { IntegrationPoint & q1 = first_pair->first; if (new_numbering.exists(q1.type, q1.ghost_type)) { UInt q1_new_el = new_numbering(q1.type, q1.ghost_type)(q1.element); AKANTU_DEBUG_ASSERT(q1_new_el != UInt(-1), "A local quadrature_point " "as been removed instead of " "just being renumbered"); q1.element = q1_new_el; nb_quad = fem.getNbIntegrationPoints(q1.type, q1.ghost_type); q1.global_num = nb_quad * q1.element + q1.num_point; } IntegrationPoint & q2 = first_pair->second; if (new_numbering.exists(q2.type, q2.ghost_type)) { UInt q2_new_el = new_numbering(q2.type, q2.ghost_type)(q2.element); AKANTU_DEBUG_ASSERT(q2_new_el != UInt(-1), "A local quadrature_point " "as been removed instead of " "just being renumbered"); q2.element = q2_new_el; nb_quad = fem.getNbIntegrationPoints(q2.type, q2.ghost_type); q2.global_num = nb_quad * q2.element + q2.num_point; } } } this->grid_synchronizer->onElementsRemoved(element_list, new_numbering, event); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index c1b6e8a60..62d7ae8be 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1628 +1,1630 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Nov 24 2015 * * @brief Implementation of the common part of the material class * * @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 "material.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), name(""), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigengradu("eigen_grad_u", *this), gradu("grad_u", *this), green_strain("green_strain", *this), piola_kirchhoff_2("piola_kirchhoff_2", *this), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material - element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension); - // model.getMesh().initElementTypeMapArray(element_filter, 1, spatial_dimension, + element_filter.initialize(model.getMesh(), + _spatial_dimension = spatial_dimension); + // model.getMesh().initElementTypeMapArray(element_filter, 1, + // spatial_dimension, // false, _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), fem(model.getFEEngine()), finite_deformation(false), name(""), model(model), spatial_dimension(dim), element_filter("element_filter", id, this->memory_id), stress("stress", *this, dim, fe_engine, this->element_filter), eigengradu("eigen_grad_u", *this, dim, fe_engine, this->element_filter), gradu("gradu", *this, dim, fe_engine, this->element_filter), green_strain("green_strain", *this, dim, fe_engine, this->element_filter), piola_kirchhoff_2("poila_kirchhoff_2", *this, dim, fe_engine, this->element_filter), potential_energy("potential_energy", *this, dim, fe_engine, this->element_filter), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse_coordinates", *this, dim, fe_engine, this->element_filter), interpolation_points_matrices("interpolation points matrices", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); element_filter.initialize(mesh, _spatial_dimension = spatial_dimension); // mesh.initElementTypeMapArray(element_filter, 1, spatial_dimension, false, // _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initialize() { registerParam("rho", rho, Real(0.), _pat_parsable | _pat_modifiable, "Density"); registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation", finite_deformation, false, _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false, _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigengradu.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); potential_energy.initialize(1); this->model.registerEventHandler(*this); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if (finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if (use_previous_stress) this->piola_kirchhoff_2.initializeHistory(); this->green_strain.initialize(spatial_dimension * spatial_dimension); } if (use_previous_stress) this->stress.initializeHistory(); if (use_previous_gradu) this->gradu.initializeHistory(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { if (it->second->hasHistory()) it->second->saveCurrentValues(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ // void Material::updateResidual(GhostType ghost_type) { // AKANTU_DEBUG_IN(); // computeAllStresses(ghost_type); // assembleResidual(ghost_type); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); if (!finite_deformation) { - Array & internal_force = const_cast &>(model.getInternalForce()); + Array & internal_force = + const_cast &>(model.getInternalForce()); Mesh & mesh = fem.getMesh(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element) { const Array & shapes_derivatives = fem.getShapesDerivatives(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); - UInt nb_quadrature_points = - fem.getNbIntegrationPoints(*it, ghost_type); + UInt nb_quadrature_points = fem.getNbIntegrationPoints(*it, ghost_type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by /// @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Array * sigma_dphi_dx = new Array(nb_element * nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array & stress_vect = this->stress(*it, ghost_type); Array::matrix_iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator B = shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_sigma_it = sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element); for (UInt q = 0; q < nb_element * nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it) Bt_sigma_it->mul(*sigma, *B); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by * @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Array * int_sigma_dphi_dx = new Array( nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); fem.integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, - size_of_shapes_derivatives, *it, ghost_type, - elem_filter); + size_of_shapes_derivatives, *it, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble - model.getDOFManager().assembleElementalArrayLocalArray(*int_sigma_dphi_dx, internal_force, - *it, ghost_type, 1, elem_filter); + model.getDOFManager().assembleElementalArrayLocalArray( + *int_sigma_dphi_dx, internal_force, *it, ghost_type, -1, + elem_filter); delete int_sigma_dphi_dx; } } } else { switch (spatial_dimension) { case 1: this->assembleInternalForces<1>(ghost_type); break; case 2: this->assembleInternalForces<2>(ghost_type); break; case 3: this->assembleInternalForces<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); - for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (const auto & type : + element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.getSize() == 0) continue; Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, - spatial_dimension, type, ghost_type, - elem_filter); + spatial_dimension, type, ghost_type, + elem_filter); gradu_vect -= eigengradu(type, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation, "The Cauchy stress can only be " "computed if you are working in " "finite deformation."); - for(auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: this->computeCauchyStress<1>(type, ghost_type); break; case 2: this->computeCauchyStress<2>(type, ghost_type); break; case 3: this->computeCauchyStress<3>(type, ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::matrix_iterator gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); Array::matrix_iterator piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF(grad_u, F_tensor); this->computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model.getDisplacement(); // resizeInternalArray(gradu); UInt spatial_dimension = model.getSpatialDimension(); - for(auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); /// compute @f$\nabla u@f$ - fem.gradientOnIntegrationPoints(displacement, gradu_vect, - spatial_dimension, type, ghost_type, - elem_filter); + fem.gradientOnIntegrationPoints(displacement, gradu_vect, spatial_dimension, + type, ghost_type, elem_filter); setToSteadyState(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); - for(auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { if (finite_deformation) { switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL<1>(type, ghost_type); assembleStiffnessMatrixL2<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrixNL<2>(type, ghost_type); assembleStiffnessMatrixL2<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrixNL<3>(type, ghost_type); assembleStiffnessMatrixL2<3>(type, ghost_type); break; } } } else { switch (spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(type, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(type, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(type, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.getSize()) { - const Array & shapes_derivatives = fem.getShapesDerivatives(type, - ghost_type); + const Array & shapes_derivatives = + fem.getShapesDerivatives(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, - type, ghost_type, elem_filter); + type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array( nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); - Array * shapesd_filtered = - new Array(nb_element, dim * nb_nodes_per_element, "filtered shapesd"); + Array * shapesd_filtered = new Array( + nb_element, dim * nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim * nb_nodes_per_element, dim * nb_nodes_per_element); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) { Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, - elem_filter); + elem_filter); delete bt_d_b; - model.getDOFManager().assembleElementalMatricesToMatrix("K", "displacement", *K_e, - type, ghost_type, _symmetric, - elem_filter); + model.getDOFManager().assembleElementalMatricesToMatrix( + "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); + const Array & shapes_derivatives = + fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); // Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); // gradu_vect.resize(nb_quadrature_points * nb_element); // fem.gradientOnIntegrationPoints(model.getIncrement(), gradu_vect, // dim, type, ghost_type, &elem_filter); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array(nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); Array::matrix_iterator Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); Array::matrix_iterator piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { Matrix & Bt_S_B = *Bt_S_B_it; Matrix & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.mul(B, S); Bt_S_B.mul(Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ - Array * K_e = new Array (nb_element, - bt_s_b_size * bt_s_b_size, - "K_e"); + Array * K_e = + new Array(nb_element, bt_s_b_size * bt_s_b_size, "K_e"); - fem.integrate(*bt_s_b, *K_e, - bt_s_b_size * bt_s_b_size, - type, ghost_type, - elem_filter); + fem.integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, + elem_filter); delete bt_s_b; - model.getDOFManager().assembleElementalMatricesToMatrix("K", "displacement", *K_e, - type, ghost_type, _symmetric, - elem_filter); + model.getDOFManager().assembleElementalMatricesToMatrix( + "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); + const Array & shapes_derivatives = + fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem.gradientOnIntegrationPoints(model.getDisplacement(), gradu_vect, dim, - type, ghost_type, elem_filter); + type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array( nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin( spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim * nb_nodes_per_element, dim * nb_nodes_per_element); Array::matrix_iterator grad_u_it = gradu_vect.begin(dim, dim); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { Matrix & grad_u = *grad_u_it; Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; // transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, // nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); fem.integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, - elem_filter); + elem_filter); delete bt_d_b; - model.getDOFManager().assembleElementalMatricesToMatrix("K", "displacement", *K_e, - type, ghost_type, _symmetric, - elem_filter); + model.getDOFManager().assembleElementalMatricesToMatrix( + "K", "displacement", *K_e, type, ghost_type, _symmetric, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template void Material::assembleInternalForces(GhostType ghost_type) { +template +void Material::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); Mesh & mesh = fem.getMesh(); - for(auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.getSize() == 0) continue; UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); - Array * shapesd_filtered = - new Array(nb_element, size_of_shapes_derivatives, "filtered shapesd"); + Array * shapesd_filtered = new Array( + nb_element, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); // Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); // Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; Array * bt_s = new Array(nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); Array::matrix_iterator grad_u_it = this->gradu(type, ghost_type).begin(dim, dim); Array::matrix_iterator grad_u_end = this->gradu(type, ghost_type).end(dim, dim); Array::matrix_iterator stress_it = this->piola_kirchhoff_2(type, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { Matrix & grad_u = *grad_u_it; Matrix & r_it = *bt_s_it; Matrix & S_it = *stress_it; setCauchyStressArray(S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.mul(B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array(nb_element, bt_s_size, "r_e"); fem.integrate(*bt_s, *r_e, bt_s_size, type, ghost_type, elem_filter); delete bt_s; - model.getDOFManager().assembleElementalArrayLocalArray(*r_e, internal_force, - type, ghost_type, 1, elem_filter); + model.getDOFManager().assembleElementalArrayLocalArray( + *r_e, internal_force, type, ghost_type, -1., elem_filter); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); - for(auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { switch (spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(type, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(type, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(type, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); Array & disp = model.getDisplacement(); fem.gradientOnIntegrationPoints(disp, gradu_vect, dim, type, ghost_type, - elem_filter); + elem_filter); UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); Array * tangent_moduli_tensors = new Array( nb_element * nb_quadrature_points, tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); tangent_moduli_tensors->clear(); computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); - Array * shapesd_filtered = - new Array(nb_element, dim * nb_nodes_per_element, "filtered shapesd"); + Array * shapesd_filtered = new Array( + nb_element, dim * nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(fem.getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array filtered_u(nb_element, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(fem.getMesh(), disp, filtered_u, type, ghost_type, elem_filter); /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator D_it = tangent_moduli_tensors->begin(tangent_moduli_size, tangent_moduli_size); Array::matrix_iterator sigma_it = stress(type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::vector_iterator u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); Matrix B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element); Vector Bu(tangent_moduli_size); Vector DBu(tangent_moduli_size); for (UInt e = 0; e < nb_element; ++e, ++u_it) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { Vector & u = *u_it; Matrix & sigma = *sigma_it; Matrix & D = *D_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bu.mul(B, u); DBu.mul(D, Bu); // Voigt notation to full symmetric tensor for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i); if (dim == 2) { sigma(0, 1) = sigma(1, 0) = DBu(2); } else if (dim == 3) { sigma(1, 2) = sigma(2, 1) = DBu(3); sigma(0, 2) = sigma(2, 0) = DBu(4); sigma(0, 1) = sigma(1, 0) = DBu(5); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); - for(auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { + for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { computePotentialEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType, GhostType) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements - for(auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { + for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { epot += fem.integrate(potential_energy(type, _not_ghost), type, _not_ghost, - element_filter(type, _not_ghost)); + element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(fem.getNbIntegrationPoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = fem.integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation( const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); this->fem.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, &(this->element_filter)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { this->fem.interpolateElementalFieldFromIntegrationPoints( this->stress, this->interpolation_points_matrices, this->interpolation_inverse_coordinates, result, ghost_type, &(this->element_filter)); } /* -------------------------------------------------------------------------- */ void Material::interpolateStressOnFacets( ElementTypeMapArray & result, ElementTypeMapArray & by_elem_result, const GhostType ghost_type) { interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model.getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); - for(auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { + for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_fil = element_filter(type, ghost_type); Array & by_elem_res = by_elem_result(type, ghost_type); UInt nb_element = elem_fil.getSize(); - UInt nb_element_full = - this->model.getMesh().getNbElement(type, ghost_type); + UInt nb_element_full = this->model.getMesh().getNbElement(type, ghost_type); UInt nb_interpolation_points_per_elem = by_elem_res.getSize() / nb_element_full; const Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); ElementType type_facet = Mesh::getFacetType(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; Element element_for_comparison(type, 0, ghost_type); - const Array > * element_to_facet = NULL; + const Array> * element_to_facet = NULL; GhostType current_ghost_type = _casper; Array * result_vec = NULL; Array::const_matrix_iterator result_it = by_elem_res.begin_reinterpret( stress_size, nb_interpolation_points_per_elem, nb_element_full); for (UInt el = 0; el < nb_element; ++el) { UInt global_el = elem_fil(el); element_for_comparison.element = global_el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element facet_elem = facet_to_element(global_el, f); UInt global_facet = facet_elem.element; if (facet_elem.ghost_type != current_ghost_type) { current_ghost_type = facet_elem.ghost_type; element_to_facet = &mesh_facets.getElementToSubelement( type_facet, current_ghost_type); result_vec = &result(type_facet, current_ghost_type); } bool is_second_element = (*element_to_facet)(global_facet)[0] != element_for_comparison; for (UInt q = 0; q < nb_quad_per_facet; ++q) { Vector result_local(result_vec->storage() + (global_facet * nb_quad_per_facet + q) * result_vec->getNbComponent() + is_second_element * stress_size, stress_size); const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template -const Array & Material::getArray(__attribute__((unused)) const ID & vect_id, - __attribute__((unused)) const ElementType & type, - __attribute__((unused)) const GhostType & ghost_type) const { +const Array & +Material::getArray(__attribute__((unused)) const ID & vect_id, + __attribute__((unused)) const ElementType & type, + __attribute__((unused)) const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(__attribute__((unused)) const ID & vect_id, __attribute__((unused)) const ElementType & type, - __attribute__((unused)) const GhostType & ghost_type) { + __attribute__((unused)) + const GhostType & ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template <> Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch (debug::Exception & e) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ template const InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template InternalField & Material::getInternal(__attribute__((unused)) const ID & int_id) { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { std::map *>::const_iterator it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { std::map *>::iterator it = internal_vectors_real.find(getID() + ":" + int_id); if (it == internal_vectors_real.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> const InternalField & Material::getInternal(const ID & int_id) const { std::map *>::const_iterator it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ template <> InternalField & Material::getInternal(const ID & int_id) { std::map *>::iterator it = internal_vectors_uint.find(getID() + ":" + int_id); if (it == internal_vectors_uint.end()) { AKANTU_SILENT_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model.getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for (; el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Array & mat_loc_num = model.getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if (el_begin == el_end) return; - ElementTypeMapArray material_local_new_numbering("remove mat filter elem", getID(), getMemoryID()); + ElementTypeMapArray material_local_new_numbering( + "remove mat filter elem", getID(), getMemoryID()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; ElementTypeMapArray::type_iterator it = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model.getMaterialLocalNumbering(type, ghost_type); if (!material_local_new_numbering.exists(type, ghost_type)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.getSize(); element.kind = (*el_begin).kind; Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if (std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter_tmp); } } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(const Array &, const NewElementsEvent &) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { UInt my_num = model.getInternalIndexFromID(getID()); - ElementTypeMapArray material_local_new_numbering("remove mat filter elem", getID(), getMemoryID()); + ElementTypeMapArray material_local_new_numbering( + "remove mat filter elem", getID(), getMemoryID()); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType gt = *g; ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); ElementTypeMapArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if (element_filter.exists(type, gt) && element_filter(type, gt).getSize()) { Array & elem_filter = element_filter(type, gt); Array & mat_indexes = this->model.getMaterialByElement(*it, gt); Array & mat_loc_num = this->model.getMaterialLocalNumbering(*it, gt); UInt nb_element = this->model.getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if (!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); Array & mat_renumbering = material_local_new_numbering(type, gt); const Array & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el; el.type = type; el.ghost_type = gt; el.kind = Mesh::getKind(type); for (UInt i = 0; i < elem_filter.getSize(); ++i) { el.element = elem_filter(i); if (std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT( new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter_tmp); } } } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeIntegrationPoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::onBeginningSolveStep(__attribute__((unused)) const AnalysisMethod & method) { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onEndSolveStep(__attribute__((unused)) const AnalysisMethod & method) { ElementTypeMapArray::type_iterator it = this->element_filter.firstType( _all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for (; it != end; ++it) { this->updateEnergies(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { ElementTypeMapArray::type_iterator it = this->element_filter.firstType( _all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for (; it != end; ++it) { this->updateEnergiesAfterDamage(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDump() { if (this->isFiniteDeformation()) this->computeAllCauchyStresses(_not_ghost); } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; std::string type = getID().substr(getID().find_last_of(":") + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, __attribute__((unused)) const Matrix & point, Matrix & extrapolated) { if (this->isInternal(id, element.kind)) { UInt nb_element = this->element_filter(element.type, element.ghost_type).getSize(); const ID name = this->getID() + ":" + id; UInt nb_quads = this->internal_vectors_real[name]->getFEEngine().getNbIntegrationPoints( element.type, element.ghost_type); const Array & internal = this->getArray(id, element.type, element.ghost_type); UInt nb_component = internal.getNbComponent(); Array::const_matrix_iterator internal_it = internal.begin_reinterpret(nb_component, nb_quads, nb_element); Element local_element = this->convertToLocalElement(element); /// instead of really extrapolating, here the value of the first GP /// is copied into the result vector. This works only for linear /// elements /// @todo extrapolate!!!! AKANTU_DEBUG_WARNING("This is a fix, values are not truly extrapolated"); const Matrix & values = internal_it[local_element.element]; UInt index = 0; Vector tmp(nb_component); for (UInt j = 0; j < values.cols(); ++j) { tmp = values(j); if (tmp.norm() > 0) { index = j; break; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } /* -------------------------------------------------------------------------- */ void Material::applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const GhostType ghost_type) { ElementTypeMapArray::type_iterator it = this->element_filter.firstType( _all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if (!element_filter(type, ghost_type).getSize()) continue; Array::matrix_iterator eigen_it = this->eigengradu(type, ghost_type) .begin(spatial_dimension, spatial_dimension); Array::matrix_iterator eigen_end = this->eigengradu(type, ghost_type) .end(spatial_dimension, spatial_dimension); for (; eigen_it != eigen_end; ++eigen_it) { Matrix & current_eigengradu = *eigen_it; current_eigengradu = prescribed_eigen_grad_u; } } } -} // akantu +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc index 8e9c52889..4c7096a16 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc @@ -1,211 +1,218 @@ /** * @file material_cohesive_bilinear.cc * * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Oct 15 2015 * * @brief Bilinear cohesive constitutive law * * @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 "material_cohesive_bilinear.hh" #include "solid_mechanics_model_cohesive.hh" +/* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ -template -MaterialCohesiveBilinear::MaterialCohesiveBilinear(SolidMechanicsModel & model, const ID & id) : - MaterialCohesiveLinear(model, id) { +template +MaterialCohesiveBilinear::MaterialCohesiveBilinear( + SolidMechanicsModel & model, const ID & id) + : MaterialCohesiveLinear(model, id) { AKANTU_DEBUG_IN(); this->registerParam("delta_0", delta_0, Real(0.), - _pat_parsable | _pat_readable, - "Elastic limit displacement"); + _pat_parsable | _pat_readable, + "Elastic limit displacement"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template +template void MaterialCohesiveBilinear::initMaterial() { AKANTU_DEBUG_IN(); this->sigma_c_eff.setRandomDistribution(this->sigma_c.getRandomParameter()); MaterialCohesiveLinear::initMaterial(); - this->delta_max .setDefaultValue(delta_0); + this->delta_max.setDefaultValue(delta_0); this->insertion_stress.setDefaultValue(0); - this->delta_max .reset(); + this->delta_max.reset(); this->insertion_stress.reset(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveBilinear::onElementsAdded(const Array & element_list, - const NewElementsEvent & event) { +template +void MaterialCohesiveBilinear::onElementsAdded( + const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); - MaterialCohesiveLinear::onElementsAdded(element_list, event); + MaterialCohesiveLinear::onElementsAdded(element_list, + event); bool scale_traction = false; // don't scale sigma_c if volume_s hasn't been specified by the user - if (!Math::are_float_equal(this->volume_s, 0.)) scale_traction = true; + if (!Math::are_float_equal(this->volume_s, 0.)) + scale_traction = true; Array::const_scalar_iterator el_it = element_list.begin(); Array::const_scalar_iterator el_end = element_list.end(); for (; el_it != el_end; ++el_it) { // filter not ghost cohesive elements - if (el_it->ghost_type != _not_ghost || el_it->kind != _ek_cohesive) continue; + if (el_it->ghost_type != _not_ghost || el_it->kind != _ek_cohesive) + continue; UInt index = el_it->element; ElementType type = el_it->type; UInt nb_element = this->model->getMesh().getNbElement(type); UInt nb_quad_per_element = this->fem_cohesive->getNbIntegrationPoints(type); - Array::vector_iterator sigma_c_begin - = this->sigma_c_eff(type).begin_reinterpret(nb_quad_per_element, nb_element); + Array::vector_iterator sigma_c_begin = + this->sigma_c_eff(type).begin_reinterpret(nb_quad_per_element, + nb_element); Vector sigma_c_vec = sigma_c_begin[index]; - Array::vector_iterator delta_c_begin - = this->delta_c_eff(type).begin_reinterpret(nb_quad_per_element, nb_element); + Array::vector_iterator delta_c_begin = + this->delta_c_eff(type).begin_reinterpret(nb_quad_per_element, + nb_element); Vector delta_c_vec = delta_c_begin[index]; - if (scale_traction) scaleTraction(*el_it, sigma_c_vec); + if (scale_traction) + scaleTraction(*el_it, sigma_c_vec); /** * Recompute sigma_c as * @f$ {\sigma_c}_\textup{new} = * \frac{{\sigma_c}_\textup{old} \delta_c} {\delta_c - \delta_0} @f$ */ for (UInt q = 0; q < nb_quad_per_element; ++q) { delta_c_vec(q) = 2 * this->G_c / sigma_c_vec(q); if (delta_c_vec(q) - delta_0 < Math::getTolerance()) - AKANTU_DEBUG_ERROR("delta_0 = " << delta_0 << - " must be lower than delta_c = " << delta_c_vec(q) - << ", modify your material file"); + AKANTU_DEBUG_ERROR("delta_0 = " + << delta_0 << " must be lower than delta_c = " + << delta_c_vec(q) << ", modify your material file"); sigma_c_vec(q) *= delta_c_vec(q) / (delta_c_vec(q) - delta_0); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveBilinear::scaleTraction(const Element & el, - Vector & sigma_c_vec) { +template +void MaterialCohesiveBilinear::scaleTraction( + const Element & el, Vector & sigma_c_vec) { AKANTU_DEBUG_IN(); Real base_sigma_c = this->sigma_c_eff; const Mesh & mesh_facets = this->model->getMeshFacets(); const FEEngine & fe_engine = this->model->getFEEngine(); - Array::const_vector_iterator coh_element_to_facet_begin - = mesh_facets.getSubelementToElement(el.type).begin(2); - const Vector & coh_element_to_facet - = coh_element_to_facet_begin[el.element]; + Array::const_vector_iterator coh_element_to_facet_begin = + mesh_facets.getSubelementToElement(el.type).begin(2); + const Vector & coh_element_to_facet = + coh_element_to_facet_begin[el.element]; // compute bounding volume Real volume = 0; // loop over facets for (UInt f = 0; f < 2; ++f) { const Element & facet = coh_element_to_facet(f); - const Array< std::vector > & facet_to_element - = mesh_facets.getElementToSubelement(facet.type, facet.ghost_type); + const Array> & facet_to_element = + mesh_facets.getElementToSubelement(facet.type, facet.ghost_type); const std::vector & element_list = facet_to_element(facet.element); std::vector::const_iterator elem = element_list.begin(); std::vector::const_iterator elem_end = element_list.end(); // loop over elements connected to each facet for (; elem != elem_end; ++elem) { // skip cohesive elements and dummy elements - if (*elem == ElementNull || elem->kind == _ek_cohesive) continue; + if (*elem == ElementNull || elem->kind == _ek_cohesive) + continue; // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); - volume += fe_engine.integrate(unit_vector, elem->type, - elem->element, elem->ghost_type); + volume += fe_engine.integrate(unit_vector, elem->type, elem->element, + elem->ghost_type); } } // scale sigma_c sigma_c_vec -= base_sigma_c; sigma_c_vec *= std::pow(this->volume_s / volume, 1. / this->m_s); sigma_c_vec += base_sigma_c; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveBilinear::computeTraction(const Array & normal, - ElementType el_type, - GhostType ghost_type) { +template +void MaterialCohesiveBilinear::computeTraction( + const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - MaterialCohesiveLinear::computeTraction(normal, - el_type, - ghost_type); + MaterialCohesiveLinear::computeTraction(normal, el_type, + ghost_type); // adjust damage - Array::scalar_iterator delta_c_it - = this->delta_c_eff(el_type, ghost_type).begin(); + Array::scalar_iterator delta_c_it = + this->delta_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator delta_max_it - = this->delta_max(el_type, ghost_type).begin(); + Array::scalar_iterator delta_max_it = + this->delta_max(el_type, ghost_type).begin(); - Array::scalar_iterator damage_it - = this->damage(el_type, ghost_type).begin(); + Array::scalar_iterator damage_it = + this->damage(el_type, ghost_type).begin(); - Array::scalar_iterator damage_end - = this->damage(el_type, ghost_type).end(); + Array::scalar_iterator damage_end = + this->damage(el_type, ghost_type).end(); for (; damage_it != damage_end; ++damage_it, ++delta_max_it, ++delta_c_it) { - *damage_it = std::max((*delta_max_it - delta_0) / (*delta_c_it - delta_0), Real(0.)); + *damage_it = + std::max((*delta_max_it - delta_0) / (*delta_c_it - delta_0), Real(0.)); *damage_it = std::min(*damage_it, Real(1.)); } } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_bilinear, MaterialCohesiveBilinear); - -} // akantu +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc index 9277b1bb1..a16acc590 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc @@ -1,700 +1,700 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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 "dof_synchronizer.hh" #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::MaterialCohesiveLinear( SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model, id), sigma_c_eff("sigma_c_eff", *this), delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this), opening_prec("opening_prec", *this), reduction_penalty("reduction_penalty", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable | _pat_readable, "Beta parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); this->registerParam("penalty", penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty coefficient"); this->registerParam("volume_s", volume_s, Real(0.), _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, Real(1.), _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa", kappa, Real(1.), _pat_parsable | _pat_readable, "Kappa parameter"); this->registerParam( "contact_after_breaking", contact_after_breaking, false, _pat_parsable | _pat_readable, "Activation of contact when the elements are fully damaged"); this->registerParam("max_quad_stress_insertion", max_quad_stress_insertion, false, _pat_parsable | _pat_readable, "Insertion of cohesive element when stress is high " "enough just on one quadrature point"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); this->use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); /// compute scalars beta2_kappa2 = beta * beta / kappa / kappa; beta2_kappa = beta * beta / kappa; if (Math::are_float_equal(beta, 0)) beta2_inv = 0; else beta2_inv = 1. / beta / beta; sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); opening_prec.initialize(spatial_dimension); reduction_penalty.initialize(1); if (!Math::are_float_equal(delta_c, 0.)) delta_c_eff.setDefaultValue(delta_c); else delta_c_eff.setDefaultValue(2 * G_c / sigma_c); if (model->getIsExtrinsic()) scaleInsertionTraction(); opening_prec.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::scaleInsertionTraction() { AKANTU_DEBUG_IN(); // do nothing if volume_s hasn't been specified by the user if (Math::are_float_equal(volume_s, 0.)) return; const Mesh & mesh_facets = model->getMeshFacets(); const FEEngine & fe_engine = model->getFEEngine(); const FEEngine & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); // loop over facet type Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); Real base_sigma_c = sigma_c; for (; first != last; ++first) { ElementType type_facet = *first; const Array > & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.getSize(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet Array::vector_iterator sigma_c_iterator = sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet); for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) { const std::vector & element_list = facet_to_element(f); // compute bounding volume Real volume = 0; std::vector::const_iterator elem = element_list.begin(); std::vector::const_iterator elem_end = element_list.end(); for (; elem != elem_end; ++elem) { if (*elem == ElementNull) continue; // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } // scale sigma_c *sigma_c_iterator -= base_sigma_c; *sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s); *sigma_c_iterator += base_sigma_c; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkInsertion( bool check_only) { AKANTU_DEBUG_IN(); const Mesh & mesh_facets = model->getMeshFacets(); CohesiveElementInserter & inserter = model->getElementInserter(); Real tolerance = Math::getTolerance(); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { ElementType type_facet = *it; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const Array & facets_check = inserter.getCheckFacets(type_facet); Array & f_insertion = inserter.getInsertionFacets(type_facet); Array & f_filter = facet_filter(type_facet); Array & sig_c_eff = sigma_c_eff(type_cohesive); Array & del_c = delta_c_eff(type_cohesive); Array & ins_stress = insertion_stress(type_cohesive); Array & trac_old = tractions_old(type_cohesive); Array & open_prec = opening_prec(type_cohesive); Array & red_penalty = reduction_penalty(type_cohesive); const Array & f_stress = model->getStressOnFacets(type_facet); const Array & sigma_lim = sigma_c(type_facet); Real max_ratio = 0.; UInt index_f = 0; UInt index_filter = 0; UInt nn = 0; UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive); AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, "The cohesive element and the corresponding facet do " "not have the same numbers of integration points"); UInt nb_facet = f_filter.getSize(); // if (nb_facet == 0) continue; Array::const_iterator sigma_lim_it = sigma_lim.begin(); Matrix stress_tmp(spatial_dimension, spatial_dimension); Matrix normal_traction(spatial_dimension, nb_quad_facet); Vector stress_check(nb_quad_facet); UInt sp2 = spatial_dimension * spatial_dimension; const Array & tangents = model->getTangents(type_facet); const Array & normals = model->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(type_facet); Array::const_vector_iterator normal_begin = normals.begin(spatial_dimension); Array::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent()); Array::const_matrix_iterator facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector > new_normal_traction; std::vector new_delta_c; // loop over each facet belonging to this material for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) { UInt facet = f_filter(f); // skip facets where check shouldn't be realized if (!facets_check(facet)) continue; // compute the effective norm on each quadrature point of the facet for (UInt q = 0; q < nb_quad_facet; ++q) { UInt current_quad = facet * nb_quad_facet + q; const Vector & normal = normal_begin[current_quad]; const Vector & tangent = tangent_begin[current_quad]; const Matrix & facet_stress_it = facet_stress_begin[current_quad]; // compute average stress on the current quadrature point Matrix stress_1(facet_stress_it.storage(), spatial_dimension, spatial_dimension); Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, spatial_dimension); stress_tmp.copy(stress_1); stress_tmp += stress_2; stress_tmp /= 2.; Vector normal_traction_vec(normal_traction(q)); // compute normal and effective stress stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent, normal_traction_vec); } // verify if the effective stress overcomes the threshold Real final_stress = stress_check.mean(); if (max_quad_stress_insertion) final_stress = *std::max_element( stress_check.storage(), stress_check.storage() + nb_quad_facet); if (final_stress > (*sigma_lim_it - tolerance)) { if (model->isDefaultSolverExplicit()) { f_insertion(facet) = true; if (!check_only) { // store the new cohesive material parameters for each quadrature // point for (UInt q = 0; q < nb_quad_facet; ++q) { Real new_sigma = stress_check(q); Vector normal_traction_vec(normal_traction(q)); if (spatial_dimension != 3) normal_traction_vec *= -1.; new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (Math::are_float_equal(delta_c, 0.)) new_delta = 2 * G_c / new_sigma; else new_delta = (*sigma_lim_it) / new_sigma * delta_c; new_delta_c.push_back(new_delta); } } } else { Real ratio = final_stress / (*sigma_lim_it); if (ratio > max_ratio) { ++nn; max_ratio = ratio; index_f = f; index_filter = f_filter(f); } } } } /// Insertion of only 1 cohesive element in case of implicit approach. The /// one subjected to the highest stress. if (!model->isDefaultSolverExplicit()) { StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Array abs_max(comm.getNbProc()); abs_max(comm.whoAmI()) = max_ratio; comm.allGather(abs_max); - Array::scalar_iterator it = + auto it = std::max_element(abs_max.begin(), abs_max.end()); Int pos = it - abs_max.begin(); if (pos != comm.whoAmI()) { AKANTU_DEBUG_OUT(); return; } if (nn) { f_insertion(index_filter) = true; if (!check_only) { // Array::iterator > normal_traction_it = // normal_traction.begin_reinterpret(nb_quad_facet, // spatial_dimension, nb_facet); Array::const_iterator sigma_lim_it = sigma_lim.begin(); for (UInt q = 0; q < nb_quad_cohesive; ++q) { Real new_sigma = (sigma_lim_it[index_f]); Vector normal_traction_vec(spatial_dimension, 0.0); new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (!Math::are_float_equal(delta_c, 0.)) new_delta = delta_c; else new_delta = 2 * G_c / (new_sigma); new_delta_c.push_back(new_delta); } } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.getSize(); UInt new_nb_quad_points = new_sigmas.size(); sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); ins_stress.resize(old_nb_quad_points + new_nb_quad_points); trac_old.resize(old_nb_quad_points + new_nb_quad_points); del_c.resize(old_nb_quad_points + new_nb_quad_points); open_prec.resize(old_nb_quad_points + new_nb_quad_points); red_penalty.resize(old_nb_quad_points + new_nb_quad_points); for (UInt q = 0; q < new_nb_quad_points; ++q) { sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; del_c(old_nb_quad_points + q) = new_delta_c[q]; red_penalty(old_nb_quad_points + q) = false; for (UInt dim = 0; dim < spatial_dimension; ++dim) { ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); open_prec(old_nb_quad_points + q, dim) = 0.; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators - Array::vector_iterator traction_it = + auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator opening_it = + auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// opening_prec is the opening of the previous step in the /// Newton-Raphson loop - Array::vector_iterator opening_prec_it = + auto opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator contact_traction_it = + auto contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator contact_opening_it = + auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); - Array::const_vector_iterator normal_it = + auto normal_it = normal.begin(spatial_dimension); - Array::vector_iterator traction_end = + auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); - Array::scalar_iterator sigma_c_it = + auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator delta_max_it = + auto delta_max_it = delta_max(el_type, ghost_type).begin(); - Array::scalar_iterator delta_max_prev_it = + auto delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); - Array::scalar_iterator delta_c_it = + auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); + auto damage_it = damage(el_type, ghost_type).begin(); - Array::vector_iterator insertion_stress_it = + auto insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); - Array::scalar_iterator reduction_penalty_it = + auto reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (!this->model->isDefaultSolverExplicit()) this->delta_max(el_type, ghost_type) .copy(this->delta_max.previous(el_type, ghost_type)); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++opening_prec_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++delta_max_prev_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTractionOnQuad( *traction_it, *delta_max_it, *delta_max_prev_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, *opening_it, *opening_prec_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkDeltaMax( GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set a predefined value to the parameter * delta_max_prev of the elements that have been inserted in the * last loading step for which convergence has not been * reached. This is done before reducing the loading and re-doing * the step. Otherwise, the updating of delta_max_prev would be * done with reference to the non-convergent solution. In this * function also other variables related to the contact and * friction behavior are correctly set. */ Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); /** * the variable "recompute" is set to true to activate the * procedure that reduce the penalty parameter for * compression. This procedure is set true only during the phase of * load_reduction, that has to be set in the maiin file. The * penalty parameter will be reduced only for the elements having * reduction_penalty = true. */ recompute = true; for (; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; /// define iterators - Array::scalar_iterator delta_max_it = + auto delta_max_it = delta_max(el_type, ghost_type).begin(); - Array::scalar_iterator delta_max_end = + auto delta_max_end = delta_max(el_type, ghost_type).end(); - Array::scalar_iterator delta_max_prev_it = + auto delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); - Array::scalar_iterator delta_c_it = + auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); - Array::vector_iterator opening_prec_it = + auto opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator opening_prec_prev_it = + auto opening_prec_prev_it = opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); /// loop on each quadrature point for (; delta_max_it != delta_max_end; ++delta_max_it, ++delta_max_prev_it, ++delta_c_it, ++opening_prec_it, ++opening_prec_prev_it) { if (*delta_max_prev_it == 0) /// elements inserted in the last incremental step, that did /// not converge *delta_max_it = *delta_c_it / 1000; else /// elements introduced in previous incremental steps, for /// which a correct value of delta_max_prev already exists *delta_max_it = *delta_max_prev_it; /// in case convergence is not reached, set opening_prec to the /// value referred to the previous incremental step *opening_prec_it = *opening_prec_prev_it; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::resetVariables( GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set the variables "recompute" and * "reduction_penalty" to false. It is called by solveStepCohesive * when convergence is reached. Such variables, in fact, have to be * false at the beginning of a new incremental step. */ Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); recompute = false; for (; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; - Array::scalar_iterator reduction_penalty_it = + auto reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); - Array::scalar_iterator reduction_penalty_end = + auto reduction_penalty_end = reduction_penalty(el_type, ghost_type).end(); /// loop on each quadrature point for (; reduction_penalty_it != reduction_penalty_end; ++reduction_penalty_it) { *reduction_penalty_it = false; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction( const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators - Array::matrix_iterator tangent_it = + auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); - Array::matrix_iterator tangent_end = + auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); - Array::const_vector_iterator normal_it = + auto normal_it = normal.begin(spatial_dimension); - Array::vector_iterator opening_it = + auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step - Array::scalar_iterator delta_max_it = + auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); - Array::scalar_iterator sigma_c_it = + auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator delta_c_it = + auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); - Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); + auto damage_it = damage(el_type, ghost_type).begin(); - Array::vector_iterator contact_opening_it = + auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); - Array::scalar_iterator reduction_penalty_it = + auto reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTangentTractionOnQuad( *tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_opening_it); // check if the tangential stiffness matrix is symmetric // for (UInt h = 0; h < spatial_dimension; ++h){ // for (UInt l = h; l < spatial_dimension; ++l){ // if (l > h){ // Real k_ls = (*tangent_it)[spatial_dimension*h+l]; // Real k_us = (*tangent_it)[spatial_dimension*l+h]; // // std::cout << "k_ls = " << k_ls << std::endl; // // std::cout << "k_us = " << k_us << std::endl; // if (std::abs(k_ls) > 1e-13 && std::abs(k_us) > 1e-13){ // Real error = std::abs((k_ls - k_us) / k_us); // if (error > 1e-10){ // std::cout << "non symmetric cohesive matrix" << std::endl; // // std::cout << "error " << error << std::endl; // } // } // } // } // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear, MaterialCohesiveLinear); } // akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index cd6ead754..f083dc51a 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1259 +1,1261 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Implementation of the SolidMechanicsModel class * * @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.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "static_communicator.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), BoundaryCondition(), f_m2a(1.0), displacement(nullptr), previous_displacement(nullptr), displacement_increment(nullptr), mass(nullptr), velocity(nullptr), acceleration(nullptr), external_force(nullptr), internal_force(nullptr), blocked_dofs(nullptr), current_position(nullptr), mass_matrix(nullptr), velocity_damping_matrix(nullptr), stiffness_matrix(nullptr), jacobian_matrix(nullptr), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), increment_flag(false), are_materials_instantiated(false) { //, pbc_synch(nullptr) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); this->mesh.registerEventHandler(*this, _ehp_solid_mechanics_model); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_material_id); this->registerSynchronizer(synchronizer, _gst_smm_mass); this->registerSynchronizer(synchronizer, _gst_smm_stress); this->registerSynchronizer(synchronizer, _gst_smm_boundary); this->registerSynchronizer(synchronizer, _gst_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); if (is_default_material_selector) { delete material_selector; material_selector = nullptr; } for (auto & internal : this->registered_internals) { delete internal.second; } // delete pbc_synch; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); this->method = smm_options.analysis_method; // initialize the vectors this->initArrays(); if (!this->hasDefaultSolver()) this->initNewSolver(this->method); // initialize pbc if (this->pbc_pair.size() != 0) this->initPBC(); // initialize the materials if (this->parser->getLastParsedFile() != "") { this->instantiateMaterials(); } //if (!smm_options.no_init_materials) { this->initMaterials(); //} // if (increment_flag) this->initBC(*this, *displacement, *displacement_increment, *external_force); // else // this->initBC(*this, *displacement, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } } return options; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; this->method = method; switch (this->method) { case _explicit_lumped_mass: { solver_name = "explicit_lumped"; tss_type = _tsst_dynamic_lumped; break; } case _explicit_consistent_mass: { solver_name = "explicit"; tss_type = _tsst_dynamic; break; } case _static: { solver_name = "static"; tss_type = _tsst_static; break; } case _implicit_dynamic: { solver_name = "implicit"; tss_type = _tsst_dynamic; break; } } if (!this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); this->setIntegrationScheme(solver_name, "displacement", options.integration_scheme_type["displacement"], options.solution_type["displacement"]); this->setDefaultSolver(solver_name); } } /* -------------------------------------------------------------------------- */ template void SolidMechanicsModel::allocNodalField(Array *& array, const ID & name) { if (array == nullptr) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":" << name; array = &(alloc(sstr_disp.str(), nb_nodes, Model::spatial_dimension, T())); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, __attribute__((unused)) NonLinearSolverType non_linear_solver_type) { DOFManager & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, "displacement"); this->allocNodalField(this->previous_displacement, "previous_displacement"); this->allocNodalField(this->displacement_increment, "displacement_increment"); this->allocNodalField(this->internal_force, "internal_force"); this->allocNodalField(this->external_force, "external_force"); this->allocNodalField(this->blocked_dofs, "blocked_dofs"); this->allocNodalField(this->current_position, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(this->velocity, "velocity"); this->allocNodalField(this->acceleration, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_static) { if (!dof_manager.hasMatrix("K")) { dof_manager.getNewMatrix("K", _symmetric); } if (!dof_manager.hasMatrix("J")) { dof_manager.getNewMatrix("J", "K"); } } } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initParallel(MeshPartition & partition, // DataAccessor * data_accessor) // { // AKANTU_DEBUG_IN(); // if (data_accessor == nullptr) // data_accessor = this; // synch_parallel = &createParallelSynch(partition, *data_accessor); // synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); // synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); // synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); // synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); // synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type, ghost_type); this->material_index.alloc(nb_element, 1, type, ghost_type); this->material_local_numbering.alloc(nb_element, 1, type, ghost_type); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initPBC() { // Model::initPBC(); // registerPBCSynchronizer(); // // as long as there are ones on the diagonal of the matrix, we can put // // boudandary true for slaves // std::map::iterator it = pbc_pair.begin(); // std::map::iterator end = pbc_pair.end(); // UInt dim = mesh.getSpatialDimension(); // while (it != end) { // for (UInt i = 0; i < dim; ++i) // (*blocked_dofs)((*it).first, i) = true; // ++it; // } // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::registerPBCSynchronizer() { // pbc_synch = new PBCSynchronizer(pbc_pair); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_uv); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_mass); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_res); // synch_registry->registerSynchronizer(*pbc_synch, _gst_for_dump); // // changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", - *this->internal_force, -1); + *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleJacobian() { this->assembleStiffnessMatrix(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as F_{int} = \int_{\Omega} N * \sigma d\Omega@f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->clear(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if(this->non_local_manager) this->non_local_manager->computeAllNonLocalStresses(); #endif // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(_gst_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(_gst_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; for (auto & material : materials) { need_to_reassemble |= material->hasStiffnessMatrixChanged(); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").clear(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) return; this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initializeUpdateResidualData() { // AKANTU_DEBUG_IN(); // UInt nb_nodes = mesh.getNbNodes(); // internal_force->resize(nb_nodes); // /// copy the forces in residual for boundary conditions // this->getDOFManager().assembleToResidual("displacement", // *this->external_force); // // start synchronization // this->asynchronousSynchronize(_gst_smm_uv); // this->waitEndSynchronize(_gst_smm_uv); // this->updateCurrentPosition(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::computeStresses() { // if (isExplicit()) { // // start synchronization // this->asynchronousSynchronize(_gst_smm_uv); // this->waitEndSynchronize(_gst_smm_uv); // // compute stresses on all local elements for each materials // std::vector::iterator mat_it; // for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { // Material & mat = **mat_it; // mat.computeAllStresses(_not_ghost); // } // #ifdef AKANTU_DAMAGE_NON_LOCAL // /* Computation of the non local part */ // this->non_local_manager->computeAllNonLocalStresses(); // #endif // } else { // std::vector::iterator mat_it; // for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { // Material & mat = **mat_it; // mat.computeAllStressesFromTangentModuli(_not_ghost); // } // } // } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // /** // * Initialize the solver and create the sparse matrices needed. // * // */ // void SolidMechanicsModel::initSolver(__attribute__((unused)) // SolverOptions & options) { // UInt nb_global_nodes = mesh.getNbGlobalNodes(); // jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", // _symmetric)); // // jacobian_matrix->buildProfile(mesh, *dof_synchronizer, // spatial_dimension); // if (!isExplicit()) { // delete stiffness_matrix; // std::stringstream sstr_sti; // sstr_sti << id << ":stiffness_matrix"; // stiffness_matrix = &(this->getDOFManager().getNewMatrix("stiffness", // _symmetric)); // } // if (solver) solver->initialize(options); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::initJacobianMatrix() { // // @todo make it more flexible: this is an ugly patch to treat the case of // non // // fix profile of the K matrix // delete jacobian_matrix; // jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", // "stiffness")); // std::stringstream sstr_solv; sstr_solv << id << ":solver"; // delete solver; // solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); // if(solver) // solver->initialize(_solver_no_options); // } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ // void SolidMechanicsModel::initImplicit(bool dynamic) { // AKANTU_DEBUG_IN(); // method = dynamic ? _implicit_dynamic : _static; // if (!increment) // setIncrementFlagOn(); // initSolver(); // // if(method == _implicit_dynamic) { // // if(integrator) delete integrator; // // integrator = new TrapezoidalRule2(); // // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { // return this->getDOFManager().getNewMatrix("velocity_damping", "jacobian"); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::implicitPred() { // AKANTU_DEBUG_IN(); // if(previous_displacement) previous_displacement->copy(*displacement); // if(method == _implicit_dynamic) // integrator->integrationSchemePred(time_step, // *displacement, // *velocity, // *acceleration, // *blocked_dofs); // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::implicitCorr() { // AKANTU_DEBUG_IN(); // if(method == _implicit_dynamic) { // integrator->integrationSchemeCorrDispl(time_step, // *displacement, // *velocity, // *acceleration, // *blocked_dofs, // *increment); // } else { // UInt nb_nodes = displacement->getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; // Real * incr_val = increment->storage(); // Real * disp_val = displacement->storage(); // bool * boun_val = blocked_dofs->storage(); // for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, // ++boun_val){ // *incr_val *= (1. - *boun_val); // *disp_val += *incr_val; // } // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::updateIncrement() { // AKANTU_DEBUG_IN(); // auto incr_it = this->displacement_increment->begin(spatial_dimension); // auto incr_end = this->displacement_increment->end(spatial_dimension); // auto disp_it = this->displacement->begin(spatial_dimension); // auto prev_disp_it = this->previous_displacement->begin(spatial_dimension); // for (; incr_it != incr_end; ++incr_it) { // *incr_it = *disp_it; // *incr_it -= *prev_disp_it; // } // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ void SolidMechanicsModel::updatePreviousDisplacement() { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_ASSERT( // previous_displacement, // "The previous displacement has to be initialized." // << " Are you working with Finite or Ineslactic deformations?"); // previous_displacement->copy(*displacement); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal(field_name, _ek_regular)) continue; for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::synchronizeBoundaries() { // AKANTU_DEBUG_IN(); // this->synchronize(_gst_smm_boundary); // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ void SolidMechanicsModel::synchronizeResidual() { // AKANTU_DEBUG_IN(); // this->synchronize(_gst_smm_res); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::setIncrementFlagOn() { // AKANTU_DEBUG_IN(); // if (!displacement_increment) { // this->allocNodalField(displacement_increment, spatial_dimension, // "displacement_increment"); // this->getDOFManager().registerDOFsIncrement("displacement", // *this->displacement_increment); // } // increment_flag = true; // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(min_dt, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const Vector & v = *v_it; const Vector & m = *m_it; Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) mv2 += v(i) * v(i) * m(i); } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { Array Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().getMatrix("M").matVecMul(*this->velocity, Mv); auto mv_it = Mv.begin(Model::spatial_dimension); auto mv_end = Mv.end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (; mv_it != mv_end; ++mv_it, ++v_it) { ekin += v_it->dot(*mv_it); } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } StaticCommunicator::getStaticCommunicator().allReduce(ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(Model::spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(Model::spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) work -= int_force(i) * incr_or_velo(i); else work += ext_force(i) * incr_or_velo(i); } } } StaticCommunicator::getStaticCommunicator().allReduce(work, _so_sum); if (this->method != _static) work *= this->getTimeStep(); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; for (auto & material : materials) energy += material->getEnergy(energy_id); /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(energy, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = this->mesh.getNbElement(type, ghost_type); if (!material_index.exists(type, ghost_type)) { this->material_index.alloc(nb_element, 1, type, ghost_type); this->material_local_numbering.alloc(nb_element, 1, type, ghost_type); } else { this->material_index(type, ghost_type).resize(nb_element); this->material_local_numbering(type, ghost_type).resize(nb_element); } } } ElementTypeMapArray filter("new_element_filter", this->getID(), this->getMemoryID()); for (auto & elem : element_list) { if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); for (auto & material : materials) material->onElementsAdded(element_list, event); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( __attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); - if (displacement) + if (displacement) { displacement->resize(nb_nodes, 0.); + ++displacement_release; + } if (mass) mass->resize(nb_nodes, 0.); if (velocity) velocity->resize(nb_nodes, 0.); if (acceleration) acceleration->resize(nb_nodes, 0.); if (external_force) external_force->resize(nb_nodes, 0.); if (internal_force) internal_force->resize(nb_nodes, 0.); if (blocked_dofs) blocked_dofs->resize(nb_nodes, 0.); if (current_position) current_position->resize(nb_nodes, 0.); if (previous_displacement) previous_displacement->resize(nb_nodes, 0.); if (displacement_increment) displacement_increment->resize(nb_nodes, 0.); for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } if(this->getDOFManager().hasMatrix("M")) { this->assembleMass(nodes_list); } if(this->getDOFManager().hasLumpedMatrix("M")) { this->assembleMassLumped(nodes_list); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) mesh.removeNodesFromArray(*mass, new_numbering); if (velocity) mesh.removeNodesFromArray(*velocity, new_numbering); if (acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if (internal_force) mesh.removeNodesFromArray(*internal_force, new_numbering); if (external_force) mesh.removeNodesFromArray(*external_force, new_numbering); if (blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) mesh.removeNodesFromArray(*displacement_increment, new_numbering); if (previous_displacement) mesh.removeNodesFromArray(*previous_displacement, new_numbering); // if (method != _explicit_lumped_mass) { // this->initSolver(); // } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass->printself(stream, indent + 2); velocity->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; for (auto & material : materials) { material->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( const GhostType & ghost_type) { for (auto & mat : materials) { try { ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id, this->memory_id); quadrature_points_coordinates.initialize( this->getFEEngine(), _spatial_dimension = Model::spatial_dimension, _nb_component = Model::spatial_dimension, _ghost_type = ghost_type); for (auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } auto & mat_non_local = dynamic_cast(*mat); mat_non_local.insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } catch (std::bad_cast &) { } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses( const GhostType & ghost_type) { for (auto & mat : materials) { try { auto & mat_non_local = dynamic_cast(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } catch (std::bad_cast &) { } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal(field_name, kind)) material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { try { auto & mat_non_local = dynamic_cast(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } catch (std::bad_cast &) { } } } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 49d6f509e..4da53ad32 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,815 +1,821 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Model of Solid Mechanics * * @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 "boundary_condition.hh" #include "data_accessor.hh" #include "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #include "integrator_gauss.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { struct SolidMechanicsModelOptions : public ModelOptions { explicit SolidMechanicsModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass); template SolidMechanicsModelOptions(use_named_args_t, pack &&... _pack); AnalysisMethod analysis_method; //bool no_init_materials; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public MeshEventHandler, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void initFull(named_argument::param_t && first, pack &&... _pack) { this->initFull(SolidMechanicsModelOptions{use_named_args, first, _pack...}); } /// initialize completely the model void initFull( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize the fem object needed for boundary conditions void initFEEngineBoundary(); /// register the tags associated with the parallel synchronizer // virtual void initParallel(MeshPartition * partition, // DataAccessor * data_accessor = NULL); /// allocate all vectors virtual void initArrays(); /// allocate all vectors // void initArraysPreviousDisplacment(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// init PBC synchronizer // void initPBC(); /// initialize a new solver and sets it as the default one to use void initNewSolver(const AnalysisMethod & method); /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; protected: /// allocate an array if needed template void allocNodalField(Array *& array, const ID & name); /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC // void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output // void synchronizeResidual(); protected: /// register PBC synchronizer // void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this assembles the stiffness matrix void assembleJacobian() override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ // public: // /// initialize the stuff for the explicit scheme // void initExplicit(AnalysisMethod analysis_method = // _explicit_lumped_mass); public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } // /// initialize the array needed by updateResidual (residual, // current_position) // void initializeUpdateResidualData(); protected: /// update the current position vector void updateCurrentPosition(); // /// assemble the residual for the explicit scheme // virtual void updateResidual(bool need_initialize = true); // /** // * \brief compute the acceleration from the residual // * this function is the explicit equivalent to solveDynamic in implicit // * In the case of lumped mass just divide the residual by the mass // * In the case of not lumped mass call // solveDynamic<_acceleration_corrector> // */ // void updateAcceleration(); /// Update the increment of displacement // void updateIncrement(); // /// Copy the actuel displacement into previous displacement // void updatePreviousDisplacement(); // /// Save stress and strain through EventManager // void saveStressAndStrainBeforeDamage(); // /// Update energies through EventManager // void updateEnergiesAfterDamage(); // /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix // void solveLumped(Array & x, const Array & A, // const Array & b, const Array & blocked_dofs, // Real alpha); // /// explicit integration predictor // void explicitPred(); // /// explicit integration corrector // void explicitCorr(); // public: // void solveStep(); // /* // ------------------------------------------------------------------------ // */ // /* Implicit */ // /* // ------------------------------------------------------------------------ // */ // public: // /// initialize the solver and the jacobian_matrix (called by // initImplicit) // void initSolver(); // /// initialize the stuff for the implicit solver // void initImplicit(bool dynamic = false); // /// solve Ma = f to get the initial acceleration // void initialAcceleration(); // /// assemble the stiffness matrix // void assembleStiffnessMatrix(); // public: // /** // * solve a step (predictor + convergence loop + corrector) using the // * the given convergence method (see akantu::SolveConvergenceMethod) // * and the given convergence criteria (see // * akantu::SolveConvergenceCriteria) // **/ // template // bool solveStep(Real tolerance, UInt max_iteration = 100); // template // bool solveStep(Real tolerance, Real & error, UInt max_iteration = 100, // bool do_not_factorize = false); // public: // /** // * solve Ku = f using the the given convergence method (see // * akantu::SolveConvergenceMethod) and the given convergence // * criteria (see akantu::SolveConvergenceCriteria) // **/ // template // bool solveStatic(Real tolerance, UInt max_iteration, // bool do_not_factorize = false); // /// create and return the velocity damping matrix // SparseMatrix & initVelocityDampingMatrix(); // /// implicit time integration predictor // void implicitPred(); // /// implicit time integration corrector // void implicitCorr(); // /// compute the Cauchy stress on user demand. // void computeCauchyStresses(); // // /// compute A and solve @f[ A\delta u = f_ext - f_int @f] // // template // // void solve(Array &increment, Real block_val = 1., // // bool need_factorize = true, bool has_profile_changed = // false); // protected: // /// finish the computation of residual to solve in increment // void updateResidualInternal(); // /// compute the support reaction and store it in force // void updateSupportReaction(); // private: // /// re-initialize the J matrix (to use if the profile of K changed) // void initJacobianMatrix(); /* ------------------------------------------------------------------------ */ // public: /// Update the stresses for the computation of the residual of the Stiffness /// matrix in the case of finite deformation // void computeStresses(); /// synchronize the ghost element boundaries values // void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// registers all the custom materials of a given type present in the input /// file // template void registerNewCustomMaterials(const ID & mat_type); /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); // /// Use a UIntData in the mesh to specify the material to use per // element void setMaterialIDsFromIntData(const std::string & data_name); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials void assignMaterialToElements(const ElementTypeMapArray * filter = NULL); /// reinitialize dof_synchronizer and solver (either in implicit or /// explicit) when cohesive elements are inserted // void reinitializeSolver(); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); + /// assemble the lumped mass only for a given list of nodes + void assembleMassLumped(const Array & node_list); + + /// assemble the lumped mass only for a given list of nodes + void assembleMass(const Array & node_list); + /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; inline UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: inline void splitElementByMaterial(const Array & elements, Array * elements_per_mat) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the current value of the time step // AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// void setTimeStep(Real time_step); /// return the of iterations done in the last solveStep // AKANTU_GET_MACRO(NumberIter, n_iter, UInt); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (external forces) AKANTU_GET_MACRO(Force, *external_force, Array &); /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); /// get the SolidMechnicsModel::incrementAcceleration vector // AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, // Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } inline void setMaterialSelector(MaterialSelector & selector); /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ // void setIncrementFlagOn(); // /// get the stiffness matrix // AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); // /// get the global jacobian matrix of the system // AKANTU_GET_MACRO(GlobalJacobianMatrix, *jacobian_matrix, const SparseMatrix // &); // /// get the mass matrix // AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); // /// get the velocity damping matrix // AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, // SparseMatrix &); /// get the FEEngine object to integrate or interpolate on the boundary inline FEEngine & getFEEngineBoundary(const ID & name = ""); // /// get integrator // AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder // &); // /// get synchronizer // AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const ElementSynchronizer // &); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// number of iterations // UInt n_iter; /// time step // Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement; /// increment of displacement Array * displacement_increment; /// lumped mass array Array * mass; /// velocities array Array * velocity; /// accelerations array Array * acceleration; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force; /// internal forces array Array * internal_force; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs; /// array of current position used during update residual Array * current_position; UInt current_position_release{0}; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta /// t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; #ifdef SWIGPYTHON public: #endif /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; #ifdef SWIGPYTHON protected: #endif /// class defining of to choose a material MaterialSelector * material_selector; /// define if it is the default selector or not bool is_default_material_selector; // /// integration scheme of second order used // IntegrationScheme2ndOrder *integrator; /// flag defining if the increment must be computed or not bool increment_flag; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// tells if the material are instantiated bool are_materials_instantiated; typedef std::map, ElementTypeMapArray *> flatten_internal_map; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// pointer to the pbc synchronizer // PBCSynchronizer * pbc_synch; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { typedef FromHigherDim FromStress; typedef FromSameDim FromTraction; } // namespace Neumann } // namespace BC /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" #include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc index 33ba85581..adbb273a4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc @@ -1,624 +1,615 @@ /** * @file fragment_manager.cc * * @author Marco Vocialta * * @date creation: Thu Jan 23 2014 * @date last modification: Mon Dec 14 2015 * * @brief Group manager to handle fragments * * @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 "fragment_manager.hh" #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model, bool dump_data, const ID & id, const MemoryID & memory_id) : GroupManager(model.getMesh(), id, memory_id), model(model), mass_center(0, model.getSpatialDimension(), "mass_center"), mass(0, model.getSpatialDimension(), "mass"), velocity(0, model.getSpatialDimension(), "velocity"), inertia_moments(0, model.getSpatialDimension(), "inertia_moments"), principal_directions( 0, model.getSpatialDimension() * model.getSpatialDimension(), "principal_directions"), quad_coordinates("quad_coordinates", id), mass_density("mass_density", id), nb_elements_per_fragment(0, 1, "nb_elements_per_fragment"), dump_data(dump_data) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// compute quadrature points' coordinates quad_coordinates.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(quad_coordinates, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getMesh().getNodes(), quad_coordinates); /// store mass density per quadrature point mass_density.initialize(mesh, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(mass_density, 1, spatial_dimension, // _not_ghost); storeMassDensityPerIntegrationPoint(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ class CohesiveElementFilter : public GroupManager::ClusteringFilter { public: CohesiveElementFilter(const SolidMechanicsModelCohesive & model, const Real max_damage = 1.) : model(model), is_unbroken(max_damage) {} bool operator()(const Element & el) const { if (el.kind == _ek_regular) return true; const Array & mat_indexes = model.getMaterialByElement(el.type, el.ghost_type); const Array & mat_loc_num = model.getMaterialLocalNumbering(el.type, el.ghost_type); const MaterialCohesive & mat = static_cast( model.getMaterial(mat_indexes(el.element))); UInt el_index = mat_loc_num(el.element); UInt nb_quad_per_element = model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(el.type, el.ghost_type); const Array & damage_array = mat.getDamage(el.type, el.ghost_type); AKANTU_DEBUG_ASSERT(nb_quad_per_element * el_index < damage_array.getSize(), "This quadrature point is out of range"); const Real * element_damage = damage_array.storage() + nb_quad_per_element * el_index; UInt unbroken_quads = std::count_if( element_damage, element_damage + nb_quad_per_element, is_unbroken); if (unbroken_quads > 0) return true; return false; } private: struct IsUnbrokenFunctor { IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {} bool operator()(const Real & x) { return x < max_damage; } const Real max_damage; }; const SolidMechanicsModelCohesive & model; const IsUnbrokenFunctor is_unbroken; }; /* -------------------------------------------------------------------------- */ void FragmentManager::buildFragments(Real damage_limit) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) ElementSynchronizer * cohesive_synchronizer = const_cast(model.getCohesiveSynchronizer()); if (cohesive_synchronizer) { cohesive_synchronizer->computeBufferSize(model, _gst_smmc_damage); cohesive_synchronizer->asynchronousSynchronize(model, _gst_smmc_damage); cohesive_synchronizer->waitEndSynchronize(model, _gst_smmc_damage); } #endif auto & mesh_facets = const_cast(mesh.getMeshFacets()); UInt spatial_dimension = model.getSpatialDimension(); std::string fragment_prefix("fragment"); /// generate fragments global_nb_fragment = - createClusters(spatial_dimension, fragment_prefix, - CohesiveElementFilter(model, damage_limit), &mesh_facets); + createClusters(spatial_dimension, mesh_facets, fragment_prefix, + CohesiveElementFilter(model, damage_limit)); nb_fragment = getNbElementGroups(spatial_dimension); fragment_index.resize(nb_fragment); UInt * fragment_index_it = fragment_index.storage(); /// loop over fragments for (const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { /// get fragment index std::string fragment_index_string = it->first.substr(fragment_prefix.size() + 1); std::stringstream sstr(fragment_index_string.c_str()); sstr >> *fragment_index_it; AKANTU_DEBUG_ASSERT(!sstr.fail(), "fragment_index is not an integer"); } /// compute fragments' mass computeMass(); if (dump_data) { createDumpDataArray(fragment_index, "fragments", true); createDumpDataArray(mass, "fragments mass"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeMass() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// create a unit field per quadrature point, since to compute mass /// it's neccessary to integrate only density ElementTypeMapArray unit_field("unit_field", id); unit_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); // mesh.initElementTypeMapArray(unit_field, spatial_dimension, // spatial_dimension, // _not_ghost); // ElementTypeMapArray::type_iterator it = // unit_field.firstType(spatial_dimension, _not_ghost, _ek_regular); // ElementTypeMapArray::type_iterator end = // unit_field.lastType(spatial_dimension, _not_ghost, _ek_regular); // for (; it != end; ++it) { // ElementType type = *it; // Array & field_array = unit_field(type); // UInt nb_element = mesh.getNbElement(type); // UInt nb_quad_per_element = // model.getFEEngine().getNbIntegrationPoints(type); // field_array.resize(nb_element * nb_quad_per_element); // field_array.set(1.); // } integrateFieldOnFragments(unit_field, mass); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeCenterOfMass() { AKANTU_DEBUG_IN(); /// integrate position multiplied by density integrateFieldOnFragments(quad_coordinates, mass_center); /// divide it by the fragments' mass Real * mass_storage = mass.storage(); Real * mass_center_storage = mass_center.storage(); UInt total_components = mass_center.getSize() * mass_center.getNbComponent(); for (UInt i = 0; i < total_components; ++i) mass_center_storage[i] /= mass_storage[i]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeVelocity() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); /// compute velocity per quadrature point ElementTypeMapArray velocity_field("velocity_field", id); velocity_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(velocity_field, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getVelocity(), velocity_field); /// integrate on fragments integrateFieldOnFragments(velocity_field, velocity); /// divide it by the fragments' mass Real * mass_storage = mass.storage(); Real * velocity_storage = velocity.storage(); UInt total_components = velocity.getSize() * velocity.getNbComponent(); for (UInt i = 0; i < total_components; ++i) velocity_storage[i] /= mass_storage[i]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Given the distance @f$ \mathbf{r} @f$ between a quadrature point * and its center of mass, the moment of inertia is computed as \f[ * I_\mathrm{CM} = \mathrm{tr}(\mathbf{r}\mathbf{r}^\mathrm{T}) * \mathbf{I} - \mathbf{r}\mathbf{r}^\mathrm{T} \f] for more * information check Wikipedia * (http://en.wikipedia.org/wiki/Moment_of_inertia#Identities_for_a_skew-symmetric_matrix) * */ void FragmentManager::computeInertiaMoments() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); computeCenterOfMass(); /// compute local coordinates products with respect to the center of match ElementTypeMapArray moments_coords("moments_coords", id); moments_coords.initialize(model.getFEEngine(), _nb_component = spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); // mesh.initElementTypeMapArray(moments_coords, // spatial_dimension * spatial_dimension, // spatial_dimension, _not_ghost); // /// resize the by element type // ElementTypeMapArray::type_iterator it = // moments_coords.firstType(spatial_dimension, _not_ghost, _ek_regular); // ElementTypeMapArray::type_iterator end = // moments_coords.lastType(spatial_dimension, _not_ghost, _ek_regular); // for (; it != end; ++it) { // ElementType type = *it; // Array & field_array = moments_coords(type); // UInt nb_element = mesh.getNbElement(type); // UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); // field_array.resize(nb_element * nb_quad_per_element); // } /// compute coordinates Array::const_vector_iterator mass_center_it = mass_center.begin(spatial_dimension); /// loop over fragments for (const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++mass_center_it) { const ElementTypeMapArray & el_list = it->second->getElements(); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); Array & moments_coords_array = moments_coords(type); const Array & quad_coordinates_array = quad_coordinates(type); const Array & el_list_array = el_list(type); Array::matrix_iterator moments_begin = moments_coords_array.begin(spatial_dimension, spatial_dimension); Array::const_vector_iterator quad_coordinates_begin = quad_coordinates_array.begin(spatial_dimension); Vector relative_coords(spatial_dimension); for (UInt el = 0; el < el_list_array.getSize(); ++el) { UInt global_el = el_list_array(el); /// loop over quadrature points for (UInt q = 0; q < nb_quad_per_element; ++q) { UInt global_q = global_el * nb_quad_per_element + q; Matrix moments_matrix = moments_begin[global_q]; const Vector & quad_coord_vector = quad_coordinates_begin[global_q]; /// to understand this read the documentation written just /// before this function relative_coords = quad_coord_vector; relative_coords -= *mass_center_it; moments_matrix.outerProduct(relative_coords, relative_coords); Real trace = moments_matrix.trace(); moments_matrix *= -1.; moments_matrix += Matrix::eye(spatial_dimension, trace); } } } } /// integrate moments Array integrated_moments(global_nb_fragment, spatial_dimension * spatial_dimension); integrateFieldOnFragments(moments_coords, integrated_moments); /// compute and store principal moments inertia_moments.resize(global_nb_fragment); principal_directions.resize(global_nb_fragment); Array::matrix_iterator integrated_moments_it = integrated_moments.begin(spatial_dimension, spatial_dimension); Array::vector_iterator inertia_moments_it = inertia_moments.begin(spatial_dimension); Array::matrix_iterator principal_directions_it = principal_directions.begin(spatial_dimension, spatial_dimension); for (UInt frag = 0; frag < global_nb_fragment; ++frag, ++integrated_moments_it, ++inertia_moments_it, ++principal_directions_it) { integrated_moments_it->eig(*inertia_moments_it, *principal_directions_it); } if (dump_data) createDumpDataArray(inertia_moments, "moments of inertia"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeAllData() { AKANTU_DEBUG_IN(); buildFragments(); computeVelocity(); computeInertiaMoments(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::storeMassDensityPerIntegrationPoint() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator end = mesh.lastType(spatial_dimension); for (; it != end; ++it) { ElementType type = *it; Array & mass_density_array = mass_density(type); UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); mass_density_array.resize(nb_element * nb_quad_per_element); const Array & mat_indexes = model.getMaterialByElement(type); Real * mass_density_it = mass_density_array.storage(); /// store mass_density for each element and quadrature point for (UInt el = 0; el < nb_element; ++el) { Material & mat = model.getMaterial(mat_indexes(el)); for (UInt q = 0; q < nb_quad_per_element; ++q, ++mass_density_it) *mass_density_it = mat.getRho(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::integrateFieldOnFragments( ElementTypeMapArray & field, Array & output) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); UInt nb_component = output.getNbComponent(); /// integration part output.resize(global_nb_fragment); output.clear(); UInt * fragment_index_it = fragment_index.storage(); Array::vector_iterator output_begin = output.begin(nb_component); /// loop over fragments for (const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { const ElementTypeMapArray & el_list = it->second->getElements(); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); const Array & density_array = mass_density(type); Array & field_array = field(type); const Array & elements = el_list(type); UInt nb_element = elements.getSize(); /// generate array to be integrated by filtering fragment's elements Array integration_array(elements.getSize() * nb_quad_per_element, nb_component); Array::matrix_iterator int_array_it = integration_array.begin_reinterpret(nb_quad_per_element, nb_component, nb_element); Array::matrix_iterator int_array_end = integration_array.end_reinterpret(nb_quad_per_element, nb_component, nb_element); Array::matrix_iterator field_array_begin = field_array.begin_reinterpret(nb_quad_per_element, nb_component, field_array.getSize() / nb_quad_per_element); Array::const_vector_iterator density_array_begin = density_array.begin_reinterpret(nb_quad_per_element, density_array.getSize() / nb_quad_per_element); for (UInt el = 0; int_array_it != int_array_end; ++int_array_it, ++el) { UInt global_el = elements(el); *int_array_it = field_array_begin[global_el]; /// multiply field by density const Vector & density_vector = density_array_begin[global_el]; for (UInt i = 0; i < nb_quad_per_element; ++i) { for (UInt j = 0; j < nb_component; ++j) { (*int_array_it)(i, j) *= density_vector(i); } } } /// integrate the field over the fragment Array integrated_array(elements.getSize(), nb_component); model.getFEEngine().integrate(integration_array, integrated_array, nb_component, type, _not_ghost, elements); /// sum over all elements and store the result Vector output_tmp(output_begin[*fragment_index_it]); output_tmp += std::accumulate(integrated_array.begin(nb_component), integrated_array.end(nb_component), Vector(nb_component)); } } /// sum output over all processors const StaticCommunicator & comm = mesh.getCommunicator(); comm.allReduce(output, _so_sum); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeNbElementsPerFragment() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); nb_elements_per_fragment.resize(global_nb_fragment); nb_elements_per_fragment.clear(); UInt * fragment_index_it = fragment_index.storage(); /// loop over fragments for (const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it, ++fragment_index_it) { const ElementTypeMapArray & el_list = it->second->getElements(); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_element = el_list(type).getSize(); nb_elements_per_fragment(*fragment_index_it) += nb_element; } } /// sum values over all processors const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_elements_per_fragment, _so_sum); if (dump_data) createDumpDataArray(nb_elements_per_fragment, "elements per fragment"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FragmentManager::createDumpDataArray(Array & data, std::string name, bool fragment_index_output) { AKANTU_DEBUG_IN(); if (data.getSize() == 0) return; - typedef typename Array::vector_iterator data_iterator; Mesh & mesh_not_const = const_cast(mesh); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_component = data.getNbComponent(); UInt * fragment_index_it = fragment_index.storage(); - data_iterator data_begin = data.begin(nb_component); + auto data_begin = data.begin(nb_component); /// loop over fragments - for (const_element_group_iterator it(element_group_begin()); + for (auto it = element_group_begin(); it != element_group_end(); ++it, ++fragment_index_it) { - const ElementGroup & fragment = *(it->second); + const auto & fragment = *(it->second); /// loop over cluster types - ElementGroup::type_iterator type_it = fragment.firstType(spatial_dimension); - ElementGroup::type_iterator type_end = fragment.lastType(spatial_dimension); - - for (; type_it != type_end; ++type_it) { - ElementType type = *type_it; - + for (auto & type : fragment.elementTypes(spatial_dimension)) { /// init mesh data - Array * mesh_data = mesh_not_const.getDataPointer( + auto & mesh_data = mesh_not_const.getDataPointer( name, type, _not_ghost, nb_component); - data_iterator mesh_data_begin = mesh_data->begin(nb_component); - - ElementGroup::const_element_iterator el_it = fragment.element_begin(type); - ElementGroup::const_element_iterator el_end = fragment.element_end(type); + auto mesh_data_begin = mesh_data.begin(nb_component); /// fill mesh data if (fragment_index_output) { - for (; el_it != el_end; ++el_it) { - Vector md_tmp(mesh_data_begin[*el_it]); + for (const auto & elem : fragment.getElements(type)) { + Vector md_tmp(mesh_data_begin[elem]); md_tmp(0) = *fragment_index_it; } } else { - for (; el_it != el_end; ++el_it) { - Vector md_tmp(mesh_data_begin[*el_it]); + for (const auto & elem : fragment.getElements(type)) { + Vector md_tmp(mesh_data_begin[elem]); md_tmp = data_begin[*fragment_index_it]; } } } } AKANTU_DEBUG_OUT(); } } // 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 1ae1e3240..ae687a1de 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,779 +1,788 @@ /** * @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_zip.hh" #include "dumpable_inline_impl.hh" #include "material_cohesive.hh" #include "shape_cohesive.hh" -#include - #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif - +/* -------------------------------------------------------------------------- */ +#include /* -------------------------------------------------------------------------- */ namespace akantu { const SolidMechanicsModelCohesiveOptions - default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false); + 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); + 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 = + // 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; 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.getSize() != 0) { type = *it; 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 = + // 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, + 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)) { + 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 - /// update shape functions - getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); - getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); - if (is_extrinsic) resizeFacetStress(); /// if (method != _explicit_lumped_mass) { /// this->initSolver(); /// } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SolidMechanicsModelCohesive::onNodesAdded( - const Array & doubled_nodes, const NewNodesEvent & event) { +void SolidMechanicsModelCohesive::onNodesAdded(const Array & new_nodes, + const NewNodesEvent & event) { AKANTU_DEBUG_IN(); - UInt nb_new_nodes = doubled_nodes.getSize(); - Array nodes_list(nb_new_nodes); + // Array nodes_list(nb_new_nodes); - for (UInt n = 0; n < nb_new_nodes; ++n) - nodes_list(n) = doubled_nodes(n, 1); + // for (UInt n = 0; n < nb_new_nodes; ++n) + // nodes_list(n) = doubled_nodes(n, 1); + SolidMechanicsModel::onNodesAdded(new_nodes, event); - SolidMechanicsModel::onNodesAdded(nodes_list, event); + UInt new_node, old_node; - for (UInt n = 0; n < nb_new_nodes; ++n) { + try { + const auto & cohesive_event = + dynamic_cast(event); + const auto & old_nodes = cohesive_event.getOldNodesList(); - UInt old_node = doubled_nodes(n, 0); - UInt new_node = doubled_nodes(n, 1); - - for (UInt dim = 0; dim < Model::spatial_dimension; ++dim) { - (*displacement)(new_node, dim) = (*displacement)(old_node, dim); - (*velocity)(new_node, dim) = (*velocity)(old_node, dim); - (*acceleration)(new_node, dim) = (*acceleration)(old_node, dim); - (*blocked_dofs)(new_node, dim) = (*blocked_dofs)(old_node, dim); + 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); + } + }; - if (current_position) - (*current_position)(new_node, dim) = (*current_position)(old_node, dim); + for (auto pair : zip(new_nodes, old_nodes)) { + std::tie(new_node, old_node) = pair; - // if (increment_acceleration) - // (*increment_acceleration)(new_node, dim) = - // (*increment_acceleration)(old_node, dim); + copy(*displacement); + copy(*velocity); + copy(*acceleration); + copy(*blocked_dofs); - // if (increment) - // (*increment)(new_node, dim) = (*increment)(old_node, dim); + if (current_position) + copy(*current_position); if (previous_displacement) - (*previous_displacement)(new_node, dim) = - (*previous_displacement)(old_node, dim); + 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; 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); + 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/test/test_model/test_non_local_toolbox/test_material.cc b/test/test_model/test_non_local_toolbox/test_material.cc index 2e3b88433..2da95be96 100644 --- a/test/test_model/test_non_local_toolbox/test_material.cc +++ b/test/test_model/test_non_local_toolbox/test_material.cc @@ -1,61 +1,61 @@ /** * @file test_material.cc * * @author Aurelia Isabel Cuba Ramos * * @date creation: Sat Sep 26 2015 * @date last modification: Wed Nov 25 2015 * * @brief Implementation of test material for the non-local neighborhood base * test * * @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 "test_material.hh" /* -------------------------------------------------------------------------- */ template TestMaterial::TestMaterial(SolidMechanicsModel & model, const ID & id) : Parent(model, id), grad_u_nl("grad_u non local", *this) { this->is_non_local = true; this->grad_u_nl.initialize(dim * dim); - this->model.registerNonLocalVariable( + this->model.getNonLocalManager().registerNonLocalVariable( this->gradu.getName(), grad_u_nl.getName(), dim * dim); } /* -------------------------------------------------------------------------- */ template void TestMaterial::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); - this->model.getNeighborhood("test_region") + this->model.getNonLocalManager().getNeighborhood("test_region") .registerNonLocalVariable(grad_u_nl.getName()); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // Instantiate the material for the 3 dimensions INSTANTIATE_MATERIAL(test_material, TestMaterial); /* -------------------------------------------------------------------------- */ diff --git a/test/test_model/test_non_local_toolbox/test_material_damage.cc b/test/test_model/test_non_local_toolbox/test_material_damage.cc index 89c7cfdd9..1bbfb2988 100644 --- a/test/test_model/test_non_local_toolbox/test_material_damage.cc +++ b/test/test_model/test_non_local_toolbox/test_material_damage.cc @@ -1,60 +1,60 @@ /** * @file test_material_damage.cc * * @author Aurelia Isabel Cuba Ramos * * @date creation: Sat Sep 26 2015 * @date last modification: Thu Oct 15 2015 * * @brief Implementation of test material damage * * @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 "test_material_damage.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template TestMaterialDamage::TestMaterialDamage(SolidMechanicsModel & model, const ID & id) : Parent(model, id), grad_u_nl("grad_u non local", *this) { this->is_non_local = true; this->grad_u_nl.initialize(dim * dim); - this->model.registerNonLocalVariable(this->gradu.getName(), + this->model.getNonLocalManager().registerNonLocalVariable(this->gradu.getName(), grad_u_nl.getName(), dim * dim); } /* -------------------------------------------------------------------------- */ template void TestMaterialDamage::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); Parent::initMaterial(); - this->model.getNeighborhood(this->name) + this->model.getNonLocalManager().getNeighborhood(this->name) .registerNonLocalVariable(grad_u_nl.getName()); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // Instantiate the material for the 3 dimensions INSTANTIATE_MATERIAL(test_material, TestMaterialDamage); /* -------------------------------------------------------------------------- */ 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 cf5ee2d17..5ab2f4031 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,127 +1,138 @@ /** * @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 /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; -int main(int argc, char *argv[]) { +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.initFull( + SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); model.limitInsertion(_y, -0.30, -0.20); model.updateAutomaticInsertion(); - Real time_step = model.getStableTimeStep()*0.05; + 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.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); + displacement(n, 1) += disp_update * position(n, 1); } model.checkCohesiveStress(); model.solveStep(); - if(s % 100 == 0) { + if (s % 100 == 0) { std::cout << "passing step " << s << "/" << max_steps << std::endl; } + model.dump(); } - 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_solver/test_sparse_solver_mumps.cc b/test/test_solver/test_sparse_solver_mumps.cc index 38888b508..b8ff41186 100644 --- a/test/test_solver/test_sparse_solver_mumps.cc +++ b/test/test_solver/test_sparse_solver_mumps.cc @@ -1,159 +1,159 @@ /** * @file test_sparse_matrix_product.cc * * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Sun Oct 19 2014 * * @brief test the matrix vector product in parallel * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_synchronizer.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "mesh_partition_scotch.hh" #include "sparse_matrix_aij.hh" #include "sparse_solver_mumps.hh" #include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); const UInt spatial_dimension = 1; const UInt nb_global_dof = 11; StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); // Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); Mesh mesh(spatial_dimension); if (prank == 0) { genMesh(mesh, nb_global_dof); RandomGenerator::seed(1496137735); } else { RandomGenerator::seed(2992275470); } mesh.distribute(); UInt node = 0; for (auto pos : mesh.getNodes()) { std::cout << prank << " " << node << " pos: " << pos << " [" << mesh.getNodeGlobalId(node) << "] " << mesh.getNodeType(node) << std::endl; ++node; } UInt nb_nodes = mesh.getNbNodes(); DOFManagerDefault dof_manager(mesh, "test_dof_manager"); Array x(nb_nodes); dof_manager.registerDOFs("x", x, _dst_nodal); - Array global_equation_number(nb_nodes); - dof_manager.getEquationsNumbers("x", global_equation_number); + const auto & local_equation_number = dof_manager.getEquationsNumbers("x"); auto & A = dof_manager.getNewMatrix("A", _symmetric); Array b(nb_nodes); TermsToAssemble terms; for (UInt i = 0; i < nb_nodes; ++i) { if (dof_manager.isLocalOrMasterDOF(i)) { - UInt gi = global_equation_number(i); + auto li = local_equation_number(i); + auto gi = dof_manager.localToGlobalEquationNumber(li); terms(i, i) = 1. / (1. + gi); } } dof_manager.assemblePreassembledMatrix("x", "x", "A", terms); std::stringstream str; str << "Matrix_" << prank << ".mtx"; A.saveMatrix(str.str()); for (UInt n = 0; n < nb_nodes; ++n) { b(n) = 1.; } SparseSolverMumps solver(dof_manager, "A"); solver.solve(x, b); auto & sync = dynamic_cast(dof_manager).getSynchronizer(); if (prank == 0) { Array x_gathered(dof_manager.getSystemSize()); sync.gather(x, x_gathered); debug::setDebugLevel(dblTest); std::cout << x_gathered << std::endl; debug::setDebugLevel(dblWarning); UInt d = 1.; for (auto x : x_gathered) { if (std::abs(x - d) / d > 1e-15) AKANTU_EXCEPTION("Error in the solution: " << x << " != " << d << " [" << (std::abs(x - d) / d) << "]."); ++d; } } else { sync.gather(x); } finalize(); return 0; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } }