diff --git a/src/fe_engine/shape_functions_inline_impl.cc b/src/fe_engine/shape_functions_inline_impl.cc index 1964c10e2..5a02a4c25 100644 --- a/src/fe_engine/shape_functions_inline_impl.cc +++ b/src/fe_engine/shape_functions_inline_impl.cc @@ -1,585 +1,586 @@ /** * @file shape_functions_inline_impl.cc * * @author Guillaume Anciaux * @author Fabian Barras * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Oct 27 2010 * @date last modification: Thu Oct 15 2015 * * @brief ShapeFunctions inline implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ __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::setIntegrationPointsByType(const Matrix & points, const GhostType & ghost_type) { integration_points(type, ghost_type).shallowCopy(points); } /* -------------------------------------------------------------------------- */ inline void ShapeFunctions::initElementalFieldInterpolationFromIntegrationPoints(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) \ initElementalFieldInterpolationFromIntegrationPoints(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::initElementalFieldInterpolationFromIntegrationPoints(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 = 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_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_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::interpolateElementalFieldFromIntegrationPoints(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) \ interpolateElementalFieldFromIntegrationPoints(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::interpolateElementalFieldFromIntegrationPoints(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 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::interpolateElementalFieldOnIntegrationPoints(const Array &u_el, Array &uq, GhostType ghost_type, const Array & shapes, const Array & filter_elements) const { UInt nb_element; UInt nb_points = integration_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::gradientElementalFieldOnIntegrationPoints(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 = integration_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/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc index 3900f1823..7ac2e8459 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc @@ -1,716 +1,714 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::MaterialCohesiveLinear(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id), sigma_c_eff("sigma_c_eff", *this), delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this), opening_prec("opening_prec", *this), reduction_penalty("reduction_penalty", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable | _pat_readable, "Beta parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); this->registerParam("penalty", penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty coefficient"); this->registerParam("volume_s", volume_s, Real(0.), _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, Real(1.), _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa", kappa, Real(1.), _pat_parsable | _pat_readable, "Kappa parameter"); this->registerParam("contact_after_breaking", contact_after_breaking, false, _pat_parsable | _pat_readable, "Activation of contact when the elements are fully damaged"); this->registerParam("max_quad_stress_insertion", max_quad_stress_insertion, false, _pat_parsable | _pat_readable, "Insertion of cohesive element when stress is high enough just on one quadrature point"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); this->use_previous_delta_max = true; - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); /// compute scalars beta2_kappa2 = beta * beta/kappa/kappa; beta2_kappa = beta * beta/kappa; if (Math::are_float_equal(beta, 0)) beta2_inv = 0; else beta2_inv = 1./beta/beta; sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); opening_prec.initialize(spatial_dimension); reduction_penalty.initialize(1); if (!Math::are_float_equal(delta_c, 0.)) delta_c_eff.setDefaultValue(delta_c); else delta_c_eff.setDefaultValue(2 * G_c / sigma_c); if (model->getIsExtrinsic()) scaleInsertionTraction(); opening_prec.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::scaleInsertionTraction() { AKANTU_DEBUG_IN(); // do nothing if volume_s hasn't been specified by the user if (Math::are_float_equal(volume_s, 0.)) return; const Mesh & mesh_facets = model->getMeshFacets(); const FEEngine & fe_engine = model->getFEEngine(); const FEEngine & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); // loop over facet type Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); Real base_sigma_c = sigma_c; for(;first != last; ++first) { ElementType type_facet = *first; const Array< std::vector > & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.getSize(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet Array::vector_iterator sigma_c_iterator = sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet); for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) { const std::vector & element_list = facet_to_element(f); // compute bounding volume Real volume = 0; std::vector::const_iterator elem = element_list.begin(); std::vector::const_iterator elem_end = element_list.end(); for (; elem != elem_end; ++elem) { if (*elem == ElementNull) continue; // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } // scale sigma_c *sigma_c_iterator -= base_sigma_c; *sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s); *sigma_c_iterator += base_sigma_c; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkInsertion(bool check_only) { AKANTU_DEBUG_IN(); const Mesh & mesh_facets = model->getMeshFacets(); CohesiveElementInserter & inserter = model->getElementInserter(); Real tolerance = Math::getTolerance(); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { ElementType type_facet = *it; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const Array & facets_check = inserter.getCheckFacets(type_facet); Array & f_insertion = inserter.getInsertionFacets(type_facet); Array & f_filter = facet_filter(type_facet); Array & sig_c_eff = sigma_c_eff(type_cohesive); Array & del_c = delta_c_eff(type_cohesive); Array & ins_stress = insertion_stress(type_cohesive); Array & trac_old = tractions_old(type_cohesive); Array & open_prec = opening_prec(type_cohesive); Array & red_penalty = reduction_penalty(type_cohesive); const Array & f_stress = model->getStressOnFacets(type_facet); const Array & sigma_lim = sigma_c(type_facet); Real max_ratio = 0.; UInt index_f = 0; UInt index_filter = 0; UInt nn = 0; UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); + UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine").getNbIntegrationPoints(type_cohesive); UInt nb_facet = f_filter.getSize(); // if (nb_facet == 0) continue; Array::const_iterator sigma_lim_it = sigma_lim.begin(); Matrix stress_tmp(spatial_dimension, spatial_dimension); Matrix normal_traction(spatial_dimension, nb_quad_facet); Vector stress_check(nb_quad_facet); UInt sp2 = spatial_dimension * spatial_dimension; const Array & tangents = model->getTangents(type_facet); const Array & normals = model->getFEEngine("FacetsFEEngine").getNormalsOnIntegrationPoints(type_facet); Array::const_vector_iterator normal_begin = normals.begin(spatial_dimension); Array::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent()); Array::const_matrix_iterator facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector< Vector > new_normal_traction; std::vector new_delta_c; // loop over each facet belonging to this material for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) { UInt facet = f_filter(f); // skip facets where check shouldn't be realized if (!facets_check(facet)) continue; // compute the effective norm on each quadrature point of the facet for (UInt q = 0; q < nb_quad_facet; ++q) { UInt current_quad = facet * nb_quad_facet + q; const Vector & normal = normal_begin[current_quad]; const Vector & tangent = tangent_begin[current_quad]; const Matrix & facet_stress_it = facet_stress_begin[current_quad]; // compute average stress on the current quadrature point Matrix stress_1(facet_stress_it.storage(), spatial_dimension, spatial_dimension); Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, spatial_dimension); stress_tmp.copy(stress_1); stress_tmp += stress_2; stress_tmp /= 2.; Vector normal_traction_vec(normal_traction(q)); // compute normal and effective stress stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent, normal_traction_vec); } // verify if the effective stress overcomes the threshold Real final_stress = stress_check.mean(); if (max_quad_stress_insertion) final_stress = *std::max_element(stress_check.storage(), stress_check.storage() + nb_quad_facet); if (final_stress > (*sigma_lim_it - tolerance)) { if (model->isExplicit()){ f_insertion(facet) = true; if (!check_only) { // store the new cohesive material parameters for each quadrature point for (UInt q = 0; q < nb_quad_facet; ++q) { Real new_sigma = stress_check(q); Vector normal_traction_vec(normal_traction(q)); if (spatial_dimension != 3) normal_traction_vec *= -1.; new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (Math::are_float_equal(delta_c, 0.)) new_delta = 2 * G_c / new_sigma; else new_delta = (*sigma_lim_it) / new_sigma * delta_c; new_delta_c.push_back(new_delta); } } }else{ Real ratio = final_stress/(*sigma_lim_it); if (ratio > max_ratio){ ++nn; max_ratio = ratio; index_f = f; index_filter = f_filter(f); } } } } /// Insertion of only 1 cohesive element in case of implicit approach. The one subjected to the highest stress. if (!model->isExplicit()){ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Array abs_max(comm.getNbProc()); abs_max(comm.whoAmI()) = max_ratio; comm.allGather(abs_max.storage(), 1); Array::scalar_iterator it = std::max_element(abs_max.begin(), abs_max.end()); Int pos = it - abs_max.begin(); if (pos != comm.whoAmI()) { AKANTU_DEBUG_OUT(); return; } if (nn) { f_insertion(index_filter) = true; if (!check_only) { // Array::iterator > normal_traction_it = // normal_traction.begin_reinterpret(nb_quad_facet, spatial_dimension, nb_facet); Array::const_iterator sigma_lim_it = sigma_lim.begin(); - for (UInt q = 0; q < nb_quad_facet; ++q) { - - // Vector ins_s(normal_traction_it[index_f].storage() + q * spatial_dimension, - // spatial_dimension); + for (UInt q = 0; q < nb_quad_cohesive; ++q) { Real new_sigma = (sigma_lim_it[index_f]); Vector normal_traction_vec(spatial_dimension, 0.0); new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; //set delta_c in function of G_c or a given delta_c value if (!Math::are_float_equal(delta_c, 0.)) new_delta = delta_c; else new_delta = 2 * G_c / (new_sigma); new_delta_c.push_back(new_delta); } } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.getSize(); UInt new_nb_quad_points = new_sigmas.size(); sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); ins_stress.resize(old_nb_quad_points + new_nb_quad_points); trac_old.resize(old_nb_quad_points + new_nb_quad_points); del_c.resize(old_nb_quad_points + new_nb_quad_points); open_prec.resize(old_nb_quad_points + new_nb_quad_points); red_penalty.resize(old_nb_quad_points + new_nb_quad_points); for (UInt q = 0; q < new_nb_quad_points; ++q) { sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; del_c(old_nb_quad_points + q) = new_delta_c[q]; red_penalty(old_nb_quad_points + q) = false; for (UInt dim = 0; dim < spatial_dimension; ++dim) { ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); open_prec(old_nb_quad_points + q, dim) = 0.; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// opening_prec is the opening of the previous step in the /// Newton-Raphson loop Array::vector_iterator opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator traction_end = tractions(el_type, ghost_type).end(spatial_dimension); Array::scalar_iterator sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_it = delta_max(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); Array::vector_iterator insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (! this->model->isExplicit()) this->delta_max(el_type, ghost_type).copy(this->delta_max.previous(el_type, ghost_type)); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++opening_prec_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++delta_max_prev_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTractionOnQuad(*traction_it, *delta_max_it, *delta_max_prev_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, *opening_it, *opening_prec_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkDeltaMax(GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set a predefined value to the parameter * delta_max_prev of the elements that have been inserted in the * last loading step for which convergence has not been * reached. This is done before reducing the loading and re-doing * the step. Otherwise, the updating of delta_max_prev would be * done with reference to the non-convergent solution. In this * function also other variables related to the contact and * friction behavior are correctly set. */ Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); /** * the variable "recompute" is set to true to activate the * procedure that reduce the penalty parameter for * compression. This procedure is set true only during the phase of * load_reduction, that has to be set in the maiin file. The * penalty parameter will be reduced only for the elements having * reduction_penalty = true. */ recompute = true; for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; /// define iterators Array::scalar_iterator delta_max_it = delta_max(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_end = delta_max(el_type, ghost_type).end(); Array::scalar_iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::vector_iterator opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_prec_prev_it = opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); /// loop on each quadrature point for (; delta_max_it != delta_max_end; ++delta_max_it, ++delta_max_prev_it, ++delta_c_it, ++opening_prec_it, ++opening_prec_prev_it) { if (*delta_max_prev_it == 0) /// elements inserted in the last incremental step, that did /// not converge *delta_max_it = *delta_c_it / 1000; else /// elements introduced in previous incremental steps, for /// which a correct value of delta_max_prev already exists *delta_max_it = *delta_max_prev_it; /// in case convergence is not reached, set opening_prec to the /// value referred to the previous incremental step *opening_prec_it = *opening_prec_prev_it; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::resetVariables(GhostType ghost_type) { AKANTU_DEBUG_IN(); /** * This function set the variables "recompute" and * "reduction_penalty" to false. It is called by solveStepCohesive * when convergence is reached. Such variables, in fact, have to be * false at the beginning of a new incremental step. */ Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); recompute = false; for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Array::scalar_iterator reduction_penalty_end = reduction_penalty(el_type, ghost_type).end(); /// loop on each quadrature point for (; reduction_penalty_it != reduction_penalty_end; ++reduction_penalty_it) { *reduction_penalty_it = false; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::matrix_iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step Array::scalar_iterator delta_max_it = delta_max.previous(el_type, ghost_type).begin(); Array::scalar_iterator sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); Array::vector_iterator contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTangentTractionOnQuad(*tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_opening_it); // check if the tangential stiffness matrix is symmetric // for (UInt h = 0; h < spatial_dimension; ++h){ // for (UInt l = h; l < spatial_dimension; ++l){ // if (l > h){ // Real k_ls = (*tangent_it)[spatial_dimension*h+l]; // Real k_us = (*tangent_it)[spatial_dimension*l+h]; // // std::cout << "k_ls = " << k_ls << std::endl; // // std::cout << "k_us = " << k_us << std::endl; // if (std::abs(k_ls) > 1e-13 && std::abs(k_us) > 1e-13){ // Real error = std::abs((k_ls - k_us) / k_us); // if (error > 1e-10){ // std::cout << "non symmetric cohesive matrix" << std::endl; // // std::cout << "error " << error << std::endl; // } // } // } // } // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialCohesiveLinear); __END_AKANTU__