diff --git a/src/fe_engine/shape_functions_inline_impl.cc b/src/fe_engine/shape_functions_inline_impl.cc index f301abced..e5f2f2130 100644 --- a/src/fe_engine/shape_functions_inline_impl.cc +++ b/src/fe_engine/shape_functions_inline_impl.cc @@ -1,584 +1,584 @@ /** * @file shape_functions_inline_impl.cc * * @author Guillaume Anciaux * @author Marco Vocialta * @author Nicolas Richart * @author Fabian Barras * * @date creation: Fri Jul 15 2011 * @date last modification: Fri Jun 13 2014 * * @brief ShapeFunctions inline implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "fe_engine.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt ShapeFunctions::getShapeSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_size = 0; #define GET_SHAPE_SIZE(type) \ shape_size = ElementClass::getShapeSize() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SIZE);// , #undef GET_SHAPE_SIZE AKANTU_DEBUG_OUT(); return shape_size; } /* -------------------------------------------------------------------------- */ inline UInt ShapeFunctions::getShapeDerivativesSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_derivatives_size = 0; #define GET_SHAPE_DERIVATIVES_SIZE(type) \ shape_derivatives_size = ElementClass::getShapeDerivativesSize() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_SIZE);// , #undef GET_SHAPE_DERIVATIVES_SIZE AKANTU_DEBUG_OUT(); return shape_derivatives_size; } /* -------------------------------------------------------------------------- */ template void ShapeFunctions::setControlPointsByType(const Matrix & points, const GhostType & ghost_type) { control_points(type, ghost_type).shallowCopy(points); } /* -------------------------------------------------------------------------- */ inline void ShapeFunctions::initElementalFieldInterpolationFromControlPoints(const ElementTypeMapArray & interpolation_points_coordinates, ElementTypeMapArray & interpolation_points_coordinates_matrices, ElementTypeMapArray & quad_points_coordinates_inv_matrices, const ElementTypeMapArray & quadrature_points_coordinates, const ElementTypeMapArray * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator it, last; if(element_filter) { it = element_filter->firstType(spatial_dimension, ghost_type); last = element_filter->lastType(spatial_dimension, ghost_type); } else { it = mesh.firstType(spatial_dimension, ghost_type); last = mesh.lastType(spatial_dimension, ghost_type); } for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array * elem_filter; if(element_filter) elem_filter = &((*element_filter)(type, ghost_type)); else elem_filter = &(empty_filter); #define AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS(type) \ initElementalFieldInterpolationFromControlPoints(interpolation_points_coordinates(type, ghost_type), \ interpolation_points_coordinates_matrices, \ quad_points_coordinates_inv_matrices, \ quadrature_points_coordinates(type, ghost_type), \ ghost_type, \ *elem_filter) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS); #undef AKANTU_INIT_ELEMENTAL_FIELD_INTERPOLATION_FROM_C_POINTS } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::initElementalFieldInterpolationFromControlPoints(const Array & interpolation_points_coordinates, ElementTypeMapArray & interpolation_points_coordinates_matrices, ElementTypeMapArray & quad_points_coordinates_inv_matrices, const Array & quadrature_points_coordinates, GhostType & ghost_type, const Array & element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); UInt nb_element_filter; if(element_filter == empty_filter) - nb_element_filter = element_filter.getSize(); - else nb_element_filter = nb_element; + nb_element_filter = nb_element; + else nb_element_filter = element_filter.getSize(); UInt nb_quad_per_element = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_interpolation_points_per_elem = interpolation_points_coordinates.getSize() / nb_element; AKANTU_DEBUG_ASSERT(interpolation_points_coordinates.getSize() % nb_element == 0, "Number of interpolation points should be a multiple of total number of elements"); if(!quad_points_coordinates_inv_matrices.exists(type, ghost_type)) quad_points_coordinates_inv_matrices.alloc(nb_element_filter, nb_quad_per_element*nb_quad_per_element, type, ghost_type); else quad_points_coordinates_inv_matrices(type, ghost_type).resize(nb_element_filter); if(!interpolation_points_coordinates_matrices.exists(type, ghost_type)) interpolation_points_coordinates_matrices.alloc(nb_element_filter, nb_interpolation_points_per_elem * nb_quad_per_element, type, ghost_type); else interpolation_points_coordinates_matrices(type, ghost_type).resize(nb_element_filter); Array & quad_inv_mat = quad_points_coordinates_inv_matrices(type, ghost_type); Array & interp_points_mat = interpolation_points_coordinates_matrices(type, ghost_type); Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); Array::const_matrix_iterator quad_coords_it = quadrature_points_coordinates.begin_reinterpret(spatial_dimension, nb_quad_per_element, - nb_element); + nb_element_filter); Array::const_matrix_iterator points_coords_begin = interpolation_points_coordinates.begin_reinterpret(spatial_dimension, nb_interpolation_points_per_elem, nb_element); Array::matrix_iterator inv_quad_coord_it = quad_inv_mat.begin(nb_quad_per_element, nb_quad_per_element); Array::matrix_iterator int_points_mat_it = interp_points_mat.begin(nb_interpolation_points_per_elem, nb_quad_per_element); /// loop over the elements of the current material and element type - for (UInt el = 0; el < nb_element; ++el, ++inv_quad_coord_it, + for (UInt el = 0; el < nb_element_filter; ++el, ++inv_quad_coord_it, ++int_points_mat_it, ++quad_coords_it) { /// matrix containing the quadrature points coordinates const Matrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationMatrix(quad_coords, quad_coord_matrix); /// invert the interpolation matrix inv_quad_coord_matrix.inverse(quad_coord_matrix); /// matrix containing the interpolation points coordinates const Matrix & points_coords = points_coords_begin[element_filter(el)]; /// matrix to store the interpolation points coordinates /// compatible with these functions Matrix & inv_points_coord_matrix = *int_points_mat_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationMatrix(points_coords, inv_points_coord_matrix); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void ShapeFunctions::buildInterpolationMatrix(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { switch (integration_order) { case 1:{ for (UInt i = 0; i < coordinates.cols(); ++i) coordMatrix(i, 0) = 1; break; } case 2:{ UInt nb_quadrature_points = coordMatrix.cols(); for (UInt i = 0; i < coordinates.cols(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(j-1, i); } break; } default:{ AKANTU_DEBUG_TO_IMPLEMENT(); break; } } } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::buildElementalFieldInterpolationMatrix(__attribute__((unused)) const Matrix & coordinates, __attribute__((unused)) Matrix & coordMatrix, __attribute__((unused)) UInt integration_order) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_segment_2>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_segment_3>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_triangle_3>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_triangle_6>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_tetrahedron_4>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_tetrahedron_10>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { buildInterpolationMatrix(coordinates, coordMatrix, integration_order); } /** * @todo Write a more efficient interpolation for quadrangles by * dropping unnecessary quadrature points * */ /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_quadrangle_4>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { if(integration_order!=ElementClassProperty<_quadrangle_4>::minimal_integration_order){ AKANTU_DEBUG_TO_IMPLEMENT(); } else { for (UInt i = 0; i < coordinates.cols(); ++i) { Real x = coordinates(0, i); Real y = coordinates(1, i); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; } } } /* -------------------------------------------------------------------------- */ template<> inline void ShapeFunctions::buildElementalFieldInterpolationMatrix<_quadrangle_8>(const Matrix & coordinates, Matrix & coordMatrix, UInt integration_order) const { if(integration_order!=ElementClassProperty<_quadrangle_8>::minimal_integration_order){ AKANTU_DEBUG_TO_IMPLEMENT(); } else { for (UInt i = 0; i < coordinates.cols(); ++i) { UInt j = 0; Real x = coordinates(0, i); Real y = coordinates(1, i); for (UInt e = 0; e <= 2; ++e) { for (UInt n = 0; n <= 2; ++n) { coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n); ++j; } } } } } /* -------------------------------------------------------------------------- */ void ShapeFunctions::interpolateElementalFieldFromControlPoints(const ElementTypeMapArray & field, const ElementTypeMapArray & interpolation_points_coordinates_matrices, const ElementTypeMapArray & quad_points_coordinates_inv_matrices, ElementTypeMapArray & result, const GhostType ghost_type, const ElementTypeMapArray * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); Mesh::type_iterator it, last; if(element_filter) { it = element_filter->firstType(spatial_dimension, ghost_type); last = element_filter->lastType(spatial_dimension, ghost_type); } else { it = mesh.firstType(spatial_dimension, ghost_type); last = mesh.lastType(spatial_dimension, ghost_type); } for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array * elem_filter; if(element_filter) elem_filter = &((*element_filter)(type, ghost_type)); else elem_filter = &(empty_filter); #define AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS(type) \ interpolateElementalFieldFromControlPoints(field(type, ghost_type), \ interpolation_points_coordinates_matrices(type, ghost_type), \ quad_points_coordinates_inv_matrices(type, ghost_type), \ result, \ ghost_type, \ *elem_filter) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS); #undef AKANTU_INTERPOLATE_ELEMENTAL_FIELD_FROM_C_POINTS } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::interpolateElementalFieldFromControlPoints(const Array & field, const Array & interpolation_points_coordinates_matrices, const Array & quad_points_coordinates_inv_matrices, ElementTypeMapArray & result, const GhostType ghost_type, const Array & element_filter) const { AKANTU_DEBUG_IN(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); - + UInt nb_quad_per_element = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_interpolation_points_per_elem = interpolation_points_coordinates_matrices.getNbComponent() / nb_quad_per_element; - + if(!result.exists(type, ghost_type)) result.alloc(nb_element*nb_interpolation_points_per_elem, field.getNbComponent(), type, ghost_type); if(element_filter != empty_filter) nb_element = element_filter.getSize(); - + Matrix coefficients(nb_quad_per_element, field.getNbComponent()); Array & result_vec = result(type, ghost_type); Array::const_matrix_iterator field_it = field.begin_reinterpret(field.getNbComponent(), nb_quad_per_element, nb_element); Array::const_matrix_iterator interpolation_points_coordinates_it = interpolation_points_coordinates_matrices.begin(nb_interpolation_points_per_elem, nb_quad_per_element); Array::matrix_iterator result_begin = result_vec.begin_reinterpret(field.getNbComponent(), nb_interpolation_points_per_elem, result_vec.getSize() / nb_interpolation_points_per_elem); Array::const_matrix_iterator inv_quad_coord_it = quad_points_coordinates_inv_matrices.begin(nb_quad_per_element, nb_quad_per_element); - /// loop over the elements of the current material and element type + /// loop over the elements of the current filter and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** * matrix containing the inversion of the quadrature points' * coordinates */ const Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// matrix containing the points' coordinates const Matrix & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the result Matrix res(result_begin[element_filter(el)]); res.mul(coefficients, coord); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ShapeFunctions::interpolateElementalFieldOnControlPoints(const Array &u_el, Array &uq, GhostType ghost_type, const Array & shapes, const Array & filter_elements) const { UInt nb_element; UInt nb_points = control_points(type, ghost_type).cols(); UInt nb_nodes_per_element = ElementClass::getShapeSize(); UInt nb_degree_of_freedom = u_el.getNbComponent() / nb_nodes_per_element; Array::const_matrix_iterator N_it; Array::const_matrix_iterator u_it; Array::matrix_iterator inter_u_it; Array * filtered_N = NULL; if(filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_N = new Array(0, shapes.getNbComponent()); FEEngine::filterElementalData(mesh, shapes, *filtered_N, type, ghost_type, filter_elements); N_it = filtered_N->begin_reinterpret(nb_nodes_per_element, nb_points, nb_element); } else { nb_element = mesh.getNbElement(type,ghost_type); N_it = shapes.begin_reinterpret(nb_nodes_per_element, nb_points, nb_element); } uq.resize(nb_element*nb_points); u_it = u_el.begin(nb_degree_of_freedom, nb_nodes_per_element); inter_u_it = uq.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element); for (UInt el = 0; el < nb_element; ++el, ++N_it, ++u_it, ++inter_u_it) { const Matrix & u = *u_it; const Matrix & N = *N_it; Matrix & inter_u = *inter_u_it; inter_u.mul(u, N); } delete filtered_N; } /* -------------------------------------------------------------------------- */ template void ShapeFunctions::gradientElementalFieldOnControlPoints(const Array &u_el, Array &out_nablauq, GhostType ghost_type, const Array & shapes_derivatives, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_points = control_points(type, ghost_type).cols(); UInt element_dimension = ElementClass::getNaturalSpaceDimension(); UInt nb_degree_of_freedom = u_el.getNbComponent() / nb_nodes_per_element; Array::const_matrix_iterator B_it; Array::const_matrix_iterator u_it; Array::matrix_iterator nabla_u_it; UInt nb_element; Array * filtered_B = NULL; if(filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_B = new Array(0, shapes_derivatives.getNbComponent()); FEEngine::filterElementalData(mesh, shapes_derivatives, *filtered_B, type, ghost_type, filter_elements); B_it = filtered_B->begin(element_dimension, nb_nodes_per_element); } else { B_it = shapes_derivatives.begin(element_dimension, nb_nodes_per_element); nb_element = mesh.getNbElement(type, ghost_type); } out_nablauq.resize(nb_element*nb_points); u_it = u_el.begin(nb_degree_of_freedom, nb_nodes_per_element); nabla_u_it = out_nablauq.begin(nb_degree_of_freedom, element_dimension); for (UInt el = 0; el < nb_element; ++el, ++u_it) { const Matrix & u = *u_it; for (UInt q = 0; q < nb_points; ++q, ++B_it, ++nabla_u_it) { const Matrix & B = *B_it; Matrix & nabla_u = *nabla_u_it; nabla_u.mul(u, B); } } delete filtered_B; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 03d13e29b..9a07c4d07 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1632 +1,1632 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_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_vector(false), 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 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(); 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); 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::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); if(!finite_deformation){ Array & residual = const_cast &>(model->getResidual()); 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->getNbQuadraturePoints(*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); delete sigma_dphi_dx; /// assemble fem->assembleArray(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete int_sigma_dphi_dx; } } } else{ switch (spatial_dimension){ case 1: this->assembleResidual<1>(ghost_type); break; case 2: this->assembleResidual<2>(ghost_type); break; case 3: this->assembleResidual<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(); Mesh::type_iterator it = fem->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = fem->getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); if (elem_filter.getSize() == 0) continue; Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ fem->gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, spatial_dimension, *it, ghost_type, elem_filter); gradu_vect -= eigengradu(*it, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, 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."); //resizeInternalArray(stress); Mesh::type_iterator it = fem->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = fem->getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) switch (spatial_dimension){ case 1: this->computeCauchyStress<1>(*it, ghost_type); break; case 2: this->computeCauchyStress<2>(*it, ghost_type); break; case 3: this->computeCauchyStress<3>(*it, 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(); Mesh::type_iterator it = fem->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = fem->getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ fem->gradientOnQuadraturePoints(displacement, gradu_vect, spatial_dimension, *it, ghost_type, elem_filter); setToSteadyState(*it, 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(); 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) { if(finite_deformation){ switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL < 1 > (*it, ghost_type); assembleStiffnessMatrixL2 < 1 > (*it, ghost_type); break; } case 2: { assembleStiffnessMatrixNL < 2 > (*it, ghost_type); assembleStiffnessMatrixL2 < 2 > (*it, ghost_type); break; } case 3: { assembleStiffnessMatrixNL < 3 > (*it, ghost_type); assembleStiffnessMatrixL2 < 3 > (*it, ghost_type); break; } } } else { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(*it, 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()) { SparseMatrix & K = const_cast(model->getStiffnessMatrix()); 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->getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem->gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, dim, 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(0, 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); delete bt_d_b; fem->assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); 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->getNbQuadraturePoints(type, ghost_type); //gradu_vect.resize(nb_quadrature_points * nb_element); // fem->gradientOnQuadraturePoints(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< dim >(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.mul < true, false > (B, S); Bt_S_B.mul < false, false > (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"); fem->integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; fem->assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); 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->getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); fem->gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, dim, 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 < true, false > (B, D); Bt_D_B.mul < false, false > (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); delete bt_d_b; fem->assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleResidual(GhostType ghost_type){ AKANTU_DEBUG_IN(); Array & residual = const_cast &> (model->getResidual()); Mesh & mesh = fem->getMesh(); Mesh::type_iterator it = element_filter.firstType(dim, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(dim, ghost_type); for (; it != last_type; ++it) { const Array & shapes_derivatives = fem->getShapesDerivatives(*it, ghost_type); Array & elem_filter = element_filter(*it, 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(*it); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, ghost_type); 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::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(*it, ghost_type).begin(dim, dim); Array::matrix_iterator grad_u_end = this->gradu(*it, ghost_type).end(dim, dim); Array::matrix_iterator stress_it = this->piola_kirchhoff_2(*it, 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 < true, false > (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, *it, ghost_type, elem_filter); delete bt_s; fem->assembleArray(*r_e, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); 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) { switch(spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(*it, 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->getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); Array & disp = model->getDisplacement(); fem->gradientOnQuadraturePoints(disp, gradu_vect, dim, type, ghost_type, 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(0, 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(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); for(; it != last_type; ++it) { computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(!potential_energy.exists(el_type, ghost_type)) { UInt nb_element = element_filter(el_type, ghost_type).getSize(); UInt nb_quadrature_points = fem->getNbQuadraturePoints(el_type, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, el_type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements Mesh::type_iterator it = element_filter.firstType(spatial_dimension); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += fem->integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _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->getNbQuadraturePoints(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->initElementalFieldInterpolationFromControlPoints(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->interpolateElementalFieldFromControlPoints(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) { - ElementTypeMapArray by_elem_result; - interpolateStress(by_elem_result, ghost_type); UInt stress_size = this->stress.getNbComponent(); const Mesh & mesh = this->model->getMesh(); const Mesh & mesh_facets = mesh.getMeshFacets(); Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = this->element_filter.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; 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_interpolation_points_per_elem = by_elem_res.getSize() / nb_element; + 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< std::vector > * 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); + 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); - Matrix result_tmp(result_it[elem_fil(el)]); + const Matrix & result_tmp(result_it[global_el]); result_local = result_tmp(f * nb_quad_per_facet + q); } } } } } /* -------------------------------------------------------------------------- */ template const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template Array & Material::getArray(const ID & vect_id, const ElementType & type, 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 << "(" < 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(const ID & int_id) const { AKANTU_DEBUG_TO_IMPLEMENT(); return NULL; } /* -------------------------------------------------------------------------- */ template InternalField & Material::getInternal(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()); 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->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeQuadraturePoints(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(__attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { 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()); 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->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::onBeginningSolveStep(const AnalysisMethod & method) { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onEndSolveStep(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) { if(!this->potential_energy.exists(*it, _not_ghost)) { UInt nb_element = this->element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = this->fem->getNbQuadraturePoints(*it, _not_ghost); this->potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } 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; } /* -------------------------------------------------------------------------- */ inline ElementTypeMap Material::getInternalDataPerElem(const ID & id, const ElementKind & element_kind, const ID & fe_engine_id) const { std::map *>::const_iterator internal_array = internal_vectors_real.find(this->getID()+":"+id); if (internal_array == internal_vectors_real.end() || internal_array->second->getElementKind() != element_kind) AKANTU_EXCEPTION("Cannot find internal field " << id << " in material " << this->name); InternalField & internal = *internal_array->second; InternalField::type_iterator it = internal.firstType(internal.getSpatialDimension(), _not_ghost, element_kind); InternalField::type_iterator last_type = internal.lastType(internal.getSpatialDimension(), _not_ghost, element_kind); ElementTypeMap res; for(; it != last_type; ++it) { UInt nb_quadrature_points = 0; nb_quadrature_points = model->getFEEngine(fe_engine_id).getNbQuadraturePoints(*it); res(*it) = internal.getNbComponent() * nb_quadrature_points; } return res; } /* -------------------------------------------------------------------------- */ void Material::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind) const { this->flattenInternalIntern(field_id, internal_flat, this->spatial_dimension, ghost_type, element_kind); } /* -------------------------------------------------------------------------- */ void Material::flattenInternalIntern(const std::string & field_id, ElementTypeMapArray & internal_flat, UInt spatial_dimension, const GhostType ghost_type, ElementKind element_kind, const ElementTypeMapArray * element_filter, const Mesh * mesh) const { typedef ElementTypeMapArray::type_iterator iterator; if(element_filter == NULL) element_filter = &(this->element_filter); if(mesh == NULL) mesh = &(this->model->mesh); iterator tit = element_filter->firstType(spatial_dimension, ghost_type, element_kind); iterator end = element_filter->lastType(spatial_dimension, ghost_type, element_kind); for (; tit != end; ++tit) { ElementType type = *tit; try { __attribute__((unused)) const Array & src_vect = this->getArray(field_id, type, ghost_type); } catch(debug::Exception & e) { continue; } const Array & src_vect = this->getArray(field_id, type, ghost_type); const Array & filter = (*element_filter)(type, ghost_type); // total number of elements for a given type UInt nb_element = mesh->getNbElement(type,ghost_type); // number of filtered elements UInt nb_element_src = filter.getSize(); // number of quadrature points per elem UInt nb_quad_per_elem = 0; // number of data per quadrature point UInt nb_data_per_quad = src_vect.getNbComponent(); if (!internal_flat.exists(type,ghost_type)) { internal_flat.alloc(nb_element * nb_quad_per_elem, nb_data_per_quad, type, ghost_type); } if (nb_element_src == 0) continue; nb_quad_per_elem = (src_vect.getSize() / nb_element_src); // number of data per element UInt nb_data = nb_quad_per_elem * src_vect.getNbComponent(); Array & dst_vect = internal_flat(type,ghost_type); dst_vect.resize(nb_element*nb_quad_per_elem); Array::const_scalar_iterator it = filter.begin(); Array::const_scalar_iterator end = filter.end(); Array::const_vector_iterator it_src = src_vect.begin_reinterpret(nb_data, nb_element_src); Array::vector_iterator it_dst = dst_vect.begin_reinterpret(nb_data, nb_element); for (; it != end ; ++it,++it_src) { it_dst[*it] = *it_src; } } }; /* -------------------------------------------------------------------------- */ /// extrapolate internal values void Material::extrapolateInternal(const ID & id, const Element & element, 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().getNbQuadraturePoints(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!!!! 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; continue; } } for (UInt i = 0; i < extrapolated.size(); ++i) { extrapolated(i) = values(index); } } else { Matrix default_values(extrapolated.rows(), extrapolated.cols(), 0.); extrapolated = default_values; } } __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 909ec2c84..48d8c0a59 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,613 +1,614 @@ /** * @file material.hh * * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_voigthelper.hh" #include "parser.hh" #include "parsable.hh" #include "data_accessor.hh" #include "internal_field.hh" #include "random_internal_field.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler, protected SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: #if __cplusplus > 199711L Material(const Material & mat) = delete; Material & operator=(const Material & mat) = delete; #endif /// Initialize material with defaults Material(SolidMechanicsModel & model, const ID & id = ""); /// Initialize material with custom mesh & fe_engine Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor virtual ~Material(); protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Function that materials can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(__attribute__((unused)) ElementType type, __attribute__((unused)) UInt index, __attribute__((unused)) Vector & epot_on_quad_points) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void updateEnergies(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { } virtual void updateEnergiesAfterDamage(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { } /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} public: /// extrapolate internal values virtual void extrapolateInternal(const ID & id, const Element & element, const Matrix & points, Matrix & extrapolated); /// compute the p-wave speed in the material virtual Real getPushWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// get a material celerity to compute the stable time step (default: is the push wave speed) virtual Real getCelerity(const Element & element) const { return getPushWaveSpeed(element); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleResidual(GhostType ghost_type); /// Operations before and after solveStep in implicit virtual void beforeSolveStep() {} virtual void afterSolveStep() {} /// save the stress in the previous_stress if needed virtual void savePreviousState(); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// add many elements at once void addElements(const Array & elements_to_add); /// remove many element at once void removeElements(const Array & elements_to_remove); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ void interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type = _not_ghost); /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points and store the * results per facet */ void interpolateStressOnFacets(ElementTypeMapArray & result, + ElementTypeMapArray & by_elem_result, const GhostType ghost_type = _not_ghost); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ void initElementalFieldInterpolation(const ElementTypeMapArray & interpolation_points_coordinates); /* ------------------------------------------------------------------------ */ /* Common part */ /* ------------------------------------------------------------------------ */ protected: /// assemble the residual template void assembleResidual(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation template void computeCauchyStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost); template inline void computeCauchyStressOnQuad(const Matrix & F, const Matrix & S, Matrix & cauchy, const Real & C33 = 1.0 ) const; template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type); /// write the stress tensor in the Voigt notation. template inline void SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect); inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// Size of the Stress matrix for the case of finite deformation see: Bathe et al, IJNME, Vol 9, 353-386, 1975 inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const; /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, 1975 template inline void setCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix); /// compute the potential energy by element void computePotentialEnergyByElements(); /// resize the intenals arrays virtual void resizeInternals(); public: /* ------------------------------------------------------------------------ */ /* Conversion functions */ /* ------------------------------------------------------------------------ */ template static inline void gradUToF (const Matrix & grad_u, Matrix & F); static inline void rightCauchy(const Matrix & F, Matrix & C); static inline void leftCauchy (const Matrix & F, Matrix & B); template static inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); template static inline void gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon); static inline Real stressToVonMises(const Matrix & stress); protected: /// converts global element to local element inline Element convertToLocalElement(const Element & global_element) const; /// converts local element to global element inline Element convertToGlobalElement(const Element & local_element) const; /// converts global quadrature point to local quadrature point inline QuadraturePoint convertToLocalPoint(const QuadraturePoint & global_point) const; /// converts local quadrature point to global quadrature point inline QuadraturePoint convertToGlobalPoint(const QuadraturePoint & local_point) const; /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ virtual void onElementsAdded(const Array & element_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void onBeginningSolveStep(const AnalysisMethod & method); virtual void onEndSolveStep(const AnalysisMethod & method); virtual void onDamageIteration(); virtual void onDamageUpdate(); virtual void onDump(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, Memory::getID(), const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// return the potential energy for the subset of elements contained by the material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained by the material virtual Real getEnergy(std::string energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray &); AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); AKANTU_GET_MACRO(FEEngine, *fem, FEEngine &); bool isNonLocal() const { return is_non_local; } template const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; template Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); template const InternalField & getInternal(const ID & id) const; template InternalField & getInternal(const ID & id); inline bool isInternal(const ID & id, const ElementKind & element_kind) const; virtual ElementTypeMap getInternalDataPerElem(const ID & id, const ElementKind & element_kind, const ID & fe_engine_id = "") const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); template inline const T & getParam(const ID & param) const; virtual void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; protected: /// internal variation of the flatten function that is more flexible and can /// be used by inherited materials to change some behavior virtual void flattenInternalIntern(const std::string & field_id, ElementTypeMapArray & internal_flat, UInt spatial_dimension, const GhostType ghost_type, ElementKind element_kind, const ElementTypeMapArray * element_filter = NULL, const Mesh * mesh = NULL) const; protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; std::map *> internal_vectors_bool; protected: /// Link to the fem object in the model FEEngine * fem; /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ElementTypeMapArray element_filter; /// stresses arrays ordered by element types InternalField stress; /// eigengrad_u arrays ordered by element types InternalField eigengradu; /// grad_u arrays ordered by element types InternalField gradu; /// Green Lagrange strain (Finite deformation) InternalField green_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types (Finite deformation) InternalField piola_kirchhoff_2; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// tell if the material need the previous strain state bool use_previous_gradu; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; /// vector that contains the names of all the internals that need to /// be transferred when material interfaces move std::vector internals_to_transfer; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type).end(this->spatial_dimension, \ this->spatial_dimension); \ \ this->stress(el_type, \ ghost_type).resize(this->gradu(el_type, \ ghost_type).getSize()); \ \ Array::iterator< Matrix > stress_it = \ this->stress(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ \ if(this->isFiniteDeformation()){ \ this->piola_kirchhoff_2(el_type, \ ghost_type).resize(this->gradu(el_type, \ ghost_type).getSize()); \ stress_it = \ this->piola_kirchhoff_2(el_type, \ ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ } \ \ for(;gradu_it != gradu_end; ++gradu_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type).end(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator sigma_it = \ this->stress(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ \ tangent_mat.resize(this->gradu(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(this->spatial_dimension); \ Array::matrix_iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } \ /* -------------------------------------------------------------------------- */ #define INSTANTIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc index 022e59333..9f63ba0d3 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc @@ -1,743 +1,746 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Tue May 08 2012 * @date last modification: Fri Sep 05 2014 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "dumpable_inline_impl.hh" #include "material_cohesive.hh" #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false, 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_distributed_synchronizer = NULL; global_connectivity = NULL; #endif delete material_selector; material_selector = new DefaultMaterialCohesiveSelector(*this); this->registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, 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_distributed_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, synch_parallel, 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]) == NULL) && 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(); mesh_facets.initElementTypeMapArray(facet_material, 1, spatial_dimension - 1); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { element.ghost_type = *gt; Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1, *gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, *gt); for(;first != last; ++first) { element.type = *first; Array & f_material = facet_material(*first, *gt); UInt nb_element = mesh_facets.getNbElement(*first, *gt); 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); #endif initAutomaticInsertion(); } else { try { /// In case of insertion of intrinsic elements along mesh physical surfaces const ParserSection & mesh_section = *(this->parser->getSubSections(_st_mesh).first); std::string cohesive_surfaces = mesh_section.getParameter("cohesive_surfaces"); initIntrinsicCohesiveMaterials(cohesive_surfaces); } catch(...) { /// Default intrinsic insertion initIntrinsicCohesiveMaterials(cohesive_index); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials(std::string cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer); #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); } #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 SolidMechanicsModel::initMaterials(); delete material_selector; material_selector = new MeshDataMaterialCohesiveSelector(*this); inserter->insertElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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(spatial_dimension, *gt, _ek_cohesive); Mesh::type_iterator last = mesh.lastType(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); #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, 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(spatial_dimension, type_ghost); Mesh::type_iterator last = mesh.lastType(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(), 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(); inserter->insertIntrinsicElements(); 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 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); mesh_facets.initElementTypeMapArray(quad_facets, spatial_dimension, spatial_dimension - 1); getFEEngine("FacetsFEEngine").interpolateOnQuadraturePoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); mesh.initElementTypeMapArray(elements_quad_facets, spatial_dimension, spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType elem_gt = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, elem_gt); Mesh::type_iterator last = mesh.lastType(spatial_dimension, elem_gt); for (; it != last; ++it) { ElementType type = *it; 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").getNbQuadraturePoints(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 < 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::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); if (need_initialize) initializeUpdateResidualData(); // f -= fint std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { try { MaterialCohesive & mat = dynamic_cast(**mat_it); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::updateResidual(false); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { try { MaterialCohesive & mat = dynamic_cast(**mat_it); 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").computeNormalsOnControlPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = spatial_dimension * (spatial_dimension - 1); mesh_facets.initElementTypeMapArray(tangents, tangent_components, spatial_dimension - 1); 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 facet_type = *it; const Array & normals = this->getFEEngine("FacetsFEEngine").getNormalsOnQuadPoints(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 (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*materials[m]); } catch(std::bad_cast&) { /// interpolate stress on facet quadrature points positions - materials[m]->interpolateStressOnFacets(facet_stress); + materials[m]->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 synch_registry->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 (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat_cohesive = dynamic_cast(*materials[m]); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch(std::bad_cast&) { } } /* if(static and extrinsic) { check max mean stresses and change inserter.getInsertionFacets(type_facet); } */ /// communicate data among processors synch_registry->synchronize(_gst_smmc_facets); /// insert cohesive elements return inserter->insertElements(); } /* -------------------------------------------------------------------------- */ 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_distributed_synchronizer != NULL) cohesive_distributed_synchronizer->computeAllBufferSizes(*this); #endif /// update shape functions getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); if (is_extrinsic) resizeFacetStress(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & doubled_nodes, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_new_nodes = doubled_nodes.getSize(); Array nodes_list(nb_new_nodes); for (UInt n = 0; n < nb_new_nodes; ++n) nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(nodes_list, event); for (UInt n = 0; n < nb_new_nodes; ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); for (UInt dim = 0; dim < 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); if (current_position) (*current_position)(new_node, dim) = (*current_position)(old_node, dim); if (increment_acceleration) (*increment_acceleration)(new_node, dim) = (*increment_acceleration)(old_node, dim); if (increment) (*increment)(new_node, dim) = (*increment)(old_node, dim); if (previous_displacement) (*previous_displacement)(new_node, dim) = (*previous_displacement)(old_node, dim); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod & method) { AKANTU_DEBUG_IN(); /****************************************************************************** This is required because the Cauchy stress is the stress measure that is used to check the insertion of cohesive elements ******************************************************************************/ std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; 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(spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(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").getNbQuadraturePoints(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 = this->spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = this->spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump(){ this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__