diff --git a/.gitignore b/.gitignore index 2b7903a21..2e0819225 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,7 @@ build .dir-locals.el TAGS third-party/blackdynamite/ third-party/cpp-array/ +*~ +release diff --git a/cmake/AkantuTestAndExamples.cmake b/cmake/AkantuTestAndExamples.cmake index e7b0ba75a..32952b1b0 100644 --- a/cmake/AkantuTestAndExamples.cmake +++ b/cmake/AkantuTestAndExamples.cmake @@ -1,319 +1,319 @@ #=============================================================================== # @file AkantuTestAndExamples.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Oct 25 2010 # @date last modification: Tue Jun 24 2014 # # @brief macros for tests and examples # # @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 . # # @section DESCRIPTION # #=============================================================================== set(AKANTU_DIFF_SCRIPT ${AKANTU_CMAKE_DIR}/akantu_diff.sh) #=============================================================================== function(manage_test_and_example et_name desc build_all label) string(TOUPPER ${et_name} upper_name) cmake_parse_arguments(manage_test_and_example "" "PACKAGE" "" ${ARGN} ) set(_activated ON) if(manage_test_and_example_PACKAGE) list(FIND ${_project}_PACKAGE_SYSTEM_PACKAGES_ON ${manage_test_and_example_PACKAGE} _ret) if(_ret EQUAL -1) set(_activated OFF) file(RELATIVE_PATH _dir ${PROJECT_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/${et_name}) list(APPEND AKANTU_TESTS_EXCLUDE_FILES /${_dir}) set(AKANTU_TESTS_EXCLUDE_FILES ${AKANTU_TESTS_EXCLUDE_FILES} CACHE INTERNAL "") endif() endif() if(NOT EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${et_name} AND _activated) # message("Example or test ${et_name} not present") return() endif() option(AKANTU_BUILD${label}${upper_name} "${desc}") mark_as_advanced(AKANTU_BUILD_${upper_name}) if(${build_all} OR NOT _activated) set(AKANTU_BUILD${label}${upper_name}_OLD ${AKANTU_BUILD${label}${upper_name}} CACHE INTERNAL "${desc}" FORCE) set(AKANTU_BUILD${label}${upper_name} ${_activated} CACHE INTERNAL "${desc}" FORCE) else() if(DEFINED AKANTU_BUILD${label}${upper_name}_OLD) set(AKANTU_BUILD${label}${upper_name} ${AKANTU_BUILD${label}${upper_name}_OLD} CACHE BOOL "${desc}" FORCE) unset(AKANTU_BUILD${label}${upper_name}_OLD CACHE) endif(DEFINED AKANTU_BUILD${label}${upper_name}_OLD) endif() if(AKANTU_BUILD${label}${upper_name}) add_subdirectory(${et_name}) endif(AKANTU_BUILD${label}${upper_name}) endfunction() #=============================================================================== # Tests #=============================================================================== if(AKANTU_TESTS) option(AKANTU_BUILD_ALL_TESTS "Build all tests" ON) mark_as_advanced(AKANTU_BUILD_ALL_TESTS) endif(AKANTU_TESTS) #=============================================================================== # Examples #=============================================================================== if(AKANTU_EXAMPLES) option(AKANTU_BUILD_ALL_EXAMPLES "Build all examples") # mark_as_advanced(AKANTU_BUILD_ALL_EXAMPLES) endif(AKANTU_EXAMPLES) #=============================================================================== macro(register_example example_name) add_executable(${example_name} ${ARGN}) target_link_libraries(${example_name} akantu ${AKANTU_EXTERNAL_LIBRARIES}) endmacro() #=============================================================================== macro(add_example example_name desc) manage_test_and_example(${example_name} ${desc} AKANTU_BUILD_ALL_EXAMPLES _EXAMPLE_ ${ARGN}) endmacro() # ============================================================================== # this should be a macro due to the enable_testing macro(add_test_tree dir) if(AKANTU_TESTS) enable_testing() include(CTest) mark_as_advanced(BUILD_TESTING) set(AKANTU_TESTS_EXCLUDE_FILES "" CACHE INTERNAL "") set(_akantu_current_parent_test ${dir} CACHE INTERNAL "Current test folder" FORCE) set(_akantu_${dir}_tests_count 0 CACHE INTERNAL "" FORCE) string(TOUPPER ${dir} _u_dir) set(AKANTU_BUILD_${_u_dir} ON CACHE INTERNAL "${desc}" FORCE) package_get_all_test_folders(_test_dirs) foreach(_dir ${_test_dirs}) add_subdirectory(${_dir}) endforeach() else() set(AKANTU_TESTS_EXCLUDE_FILES "${CMAKE_CURRENT_BINARY_DIR}/${dir}" CACHE INTERNAL "") endif() endmacro() # ============================================================================== function(add_akantu_test dir desc) set(_my_parent_dir ${_akantu_current_parent_test}) # initialize variables set(_akantu_current_parent_test ${dir} CACHE INTERNAL "Current test folder" FORCE) set(_akantu_${dir}_tests_count 0 CACHE INTERNAL "" FORCE) # set the option for this directory string(TOUPPER ${dir} _u_dir) option(AKANTU_BUILD_${_u_dir} "${desc}") mark_as_advanced(AKANTU_BUILD_${_u_dir}) # add the sub-directory add_subdirectory(${dir}) # if no test can be activated make the option disappear set(_force_deactivate_count FALSE) if(${_akantu_${dir}_tests_count} EQUAL 0) set(_force_deactivate_count TRUE) endif() # if parent off make the option disappear set(_force_deactivate_parent FALSE) string(TOUPPER ${_my_parent_dir} _u_parent_dir) if(NOT AKANTU_BUILD_${_u_parent_dir}) set(_force_deactivate_parent TRUE) endif() if(_force_deactivate_parent OR _force_deactivate_count OR AKANTU_BUILD_ALL_TESTS) if(NOT DEFINED _AKANTU_BUILD_${_u_dir}_SAVE) set(_AKANTU_BUILD_${_u_dir}_SAVE ${AKANTU_BUILD_${_u_dir}} CACHE INTERNAL "" FORCE) endif() unset(AKANTU_BUILD_${_u_dir} CACHE) if(AKANTU_BUILD_ALL_TESTS AND NOT _force_deactivate_count) set(AKANTU_BUILD_${_u_dir} ON CACHE INTERNAL "${desc}" FORCE) else() set(AKANTU_BUILD_${_u_dir} OFF CACHE INTERNAL "${desc}" FORCE) endif() else() if(DEFINED _AKANTU_BUILD_${_u_dir}_SAVE) unset(AKANTU_BUILD_${_u_dir} CACHE) set(AKANTU_BUILD_${_u_dir} ${_AKANTU_BUILD_${_u_dir}_SAVE} CACHE BOOL "${desc}") unset(_AKANTU_BUILD_${_u_dir}_SAVE CACHE) endif() endif() # adding up to the parent count math(EXPR _tmp_parent_count "${_akantu_${dir}_tests_count} + ${_akantu_${_my_parent_dir}_tests_count}") set(_akantu_${_my_parent_dir}_tests_count ${_tmp_parent_count} CACHE INTERNAL "" FORCE) # restoring the parent current dir set(_akantu_current_parent_test ${_my_parent_dir} CACHE INTERNAL "Current test folder" FORCE) endfunction() # ============================================================================== function(register_test test_name) set(multi_variables SOURCES FILES_TO_COPY DEPENDENCIES DIRECTORIES_TO_CREATE COMPILE_OPTIONS EXTRA_FILES ) cmake_parse_arguments(_register_test "" "PACKAGE" "${multi_variables}" ${ARGN} ) if(NOT _register_test_PACKAGE) message(FATAL_ERROR "No reference package was defined for the test ${test_name} in folder ${CMAKE_CURRENT_SOURCE_DIR}") endif() package_get_name(${_register_test_PACKAGE} _pkg_name) package_is_activated(${_pkg_name} _act) if(_act) math(EXPR _tmp_parent_count "${_akantu_${_akantu_current_parent_test}_tests_count} + 1") set(_akantu_${_akantu_current_parent_test}_tests_count ${_tmp_parent_count} CACHE INTERNAL "" FORCE) string(TOUPPER ${_akantu_current_parent_test} _u_parent) if(AKANTU_BUILD_${_u_parent}) - package_get_include_directories( + package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) - package_get_external_informations( + package_get_all_external_informations( AKANTU_EXTERNAL_INCLUDE_DIR AKANTU_EXTERNAL_LIBRARIES ) # set the proper includes to build most of the tests include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_DIR} ) # Register the executable to compile add_executable(${test_name} ${_register_test_SOURCES} ${_register_test_UNPARSED_ARGUMENTS}) set_property(TARGET ${test_name} APPEND PROPERTY INCLUDE_DIRECTORIES ${AKANTU_LIBRARY_INCLUDE_DIRS} ${AKANTU_EXTERNAL_INCLUDE_DIR}) target_link_libraries(${test_name} akantu ${AKANTU_EXTERNAL_LIBRARIES}) # add the extra compilation options if(_register_test_COMPILE_OPTIONS) set_target_properties(${test_name} PROPERTIES COMPILE_DEFINITIONS "${_register_test_COMPILE_OPTIONS}") endif() set(_test_all_files) # add the different dependencies (meshes, local libraries, ...) foreach(_dep ${_register_test_DEPENDENCIES}) add_dependencies(${test_name} ${_dep}) get_target_property(_dep_in_ressources ${_dep} RESSOURCES) if(_dep_in_ressources) list(APPEND _test_all_files ${_dep_in_ressources}) endif() endforeach() # copy the needed files to the build folder if(_register_test_FILES_TO_COPY) foreach(_file ${_register_test_FILES_TO_COPY}) file(COPY ${_file} DESTINATION .) list(APPEND _test_all_files ${CMAKE_CURRENT_SOURCE_DIR}/${_file}) endforeach() endif() # create the needed folders in the build folder if(_register_test_DIRECTORIES_TO_CREATE) foreach(_dir ${_register_test_DIRECTORIES_TO_CREATE}) if(IS_ABSOLUTE ${dir}) file(MAKE_DIRECTORY ${_dir}) else() file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/${_dir}) endif() endforeach() endif() # add the source files in the list of all files foreach(_file ${_register_test_SOURCES} ${_register_test_UNPARSED_ARGUMENTS} ${_register_test_EXTRA_FILES}) if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${_file}) list(APPEND _test_all_files ${CMAKE_CURRENT_SOURCE_DIR}/${_file}) else() message("The file \"${_file}\" registred by the test \"${test_name}\" does not exists") endif() endforeach() # register the test for ctest if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${test_name}.sh) file(COPY ${test_name}.sh DESTINATION .) list(APPEND _test_all_files ${CMAKE_CURRENT_SOURCE_DIR}/${test_name}.sh) add_test(${test_name}_run ${CMAKE_CURRENT_BINARY_DIR}/${test_name}.sh) elseif(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${test_name}.verified) list(APPEND _test_all_files ${CMAKE_CURRENT_SOURCE_DIR}/${test_name}.verified) add_test(${test_name}_run ${AKANTU_DIFF_SCRIPT} ${test_name} ${CMAKE_CURRENT_SOURCE_DIR}/${test_name}.verified) else() add_test(${test_name}_run ${CMAKE_CURRENT_BINARY_DIR}/${test_name}) endif() # add the executable as a dependency of the run set_tests_properties(${test_name}_run PROPERTIES DEPENDS ${test_name}) # clean the list of all files for this test and add them in the total list set(_tmp ${AKANTU_TESTS_FILES}) foreach(_file ${_source_file}) get_filename_component(_full ${_file} ABSOLUTE) file(RELATIVE_PATH __file ${PROJECT_SOURCE_DIR} ${_full}) list(APPEND _tmp ${__file}) list(APPEND _pkg_tmp ${__file}) endforeach() set(AKANTU_TESTS_FILES ${_tmp} CACHE INTERNAL "") endif() endif() endfunction() \ No newline at end of file diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc index 5927ebfca..a6a7daa84 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc @@ -1,210 +1,210 @@ /** * @file material_cohesive_bilinear.cc * * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jul 31 2014 * * @brief Bilinear cohesive constitutive law * * @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 "material_cohesive_bilinear.hh" #include "solid_mechanics_model_cohesive.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveBilinear::MaterialCohesiveBilinear(SolidMechanicsModel & model, const ID & id) : MaterialCohesiveLinear(model, id) { AKANTU_DEBUG_IN(); this->registerParam("delta_0", delta_0, 0., _pat_parsable | _pat_readable, "Elastic limit displacement"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::initMaterial() { AKANTU_DEBUG_IN(); this->sigma_c_eff.setRandomDistribution(this->sigma_c.getRandomParameter()); MaterialCohesiveLinear::initMaterial(); this->delta_max .setDefaultValue(delta_0); this->insertion_stress.setDefaultValue(0); this->delta_max .reset(); this->insertion_stress.reset(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::onElementsAdded(element_list, event); bool scale_traction = false; // don't scale sigma_c if volume_s hasn't been specified by the user if (!Math::are_float_equal(this->volume_s, 0.)) scale_traction = true; Array::const_scalar_iterator el_it = element_list.begin(); Array::const_scalar_iterator el_end = element_list.end(); for (; el_it != el_end; ++el_it) { // filter not ghost cohesive elements if (el_it->ghost_type != _not_ghost || el_it->kind != _ek_cohesive) continue; UInt index = el_it->element; ElementType type = el_it->type; UInt nb_element = this->model->getMesh().getNbElement(type); UInt nb_quad_per_element = this->fem_cohesive->getNbQuadraturePoints(type); Array::vector_iterator sigma_c_begin = this->sigma_c_eff(type).begin_reinterpret(nb_quad_per_element, nb_element); Vector & sigma_c_vec = sigma_c_begin[index]; Array::vector_iterator delta_c_begin - = this->delta_c(type).begin_reinterpret(nb_quad_per_element, nb_element); + = this->delta_c_eff(type).begin_reinterpret(nb_quad_per_element, nb_element); Vector & delta_c_vec = delta_c_begin[index]; if (scale_traction) scaleTraction(*el_it, sigma_c_vec); /** * Recompute sigma_c as * @f$ {\sigma_c}_\textup{new} = * \frac{{\sigma_c}_\textup{old} \delta_c} {\delta_c - \delta_0} @f$ */ for (UInt q = 0; q < nb_quad_per_element; ++q) { - delta_c_vec(q) = 2 * this->G_cI / sigma_c_vec(q); + delta_c_vec(q) = 2 * this->G_c / sigma_c_vec(q); if (delta_c_vec(q) - delta_0 < Math::getTolerance()) AKANTU_DEBUG_ERROR("delta_0 = " << delta_0 << " must be lower than delta_c = " << delta_c_vec(q) << ", modify your material file"); sigma_c_vec(q) *= delta_c_vec(q) / (delta_c_vec(q) - delta_0); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::scaleTraction(const Element & el, Vector & sigma_c_vec) { AKANTU_DEBUG_IN(); Real base_sigma_c = this->sigma_c_eff; const Mesh & mesh_facets = this->model->getMeshFacets(); const FEEngine & fe_engine = this->model->getFEEngine(); Array::const_vector_iterator coh_element_to_facet_begin = mesh_facets.getSubelementToElement(el.type).begin(2); const Vector & coh_element_to_facet = coh_element_to_facet_begin[el.element]; // compute bounding volume Real volume = 0; // loop over facets for (UInt f = 0; f < 2; ++f) { const Element & facet = coh_element_to_facet(f); const Array< std::vector > & facet_to_element = mesh_facets.getElementToSubelement(facet.type, facet.ghost_type); const std::vector & element_list = facet_to_element(facet.element); std::vector::const_iterator elem = element_list.begin(); std::vector::const_iterator elem_end = element_list.end(); // loop over elements connected to each facet for (; elem != elem_end; ++elem) { // skip cohesive elements and dummy elements if (*elem == ElementNull || elem->kind == _ek_cohesive) continue; // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbQuadraturePoints(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_vec -= base_sigma_c; sigma_c_vec *= std::pow(this->volume_s / volume, 1. / this->m_s); sigma_c_vec += base_sigma_c; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveBilinear::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialCohesiveLinear::computeTraction(normal, el_type, ghost_type); // adjust damage Array::scalar_iterator delta_c_it - = this->delta_c(el_type, ghost_type).begin(); + = this->delta_c_eff(el_type, ghost_type).begin(); Array::scalar_iterator delta_max_it = this->delta_max(el_type, ghost_type).begin(); Array::scalar_iterator damage_it = this->damage(el_type, ghost_type).begin(); Array::scalar_iterator damage_end = this->damage(el_type, ghost_type).end(); for (; damage_it != damage_end; ++damage_it, ++delta_max_it, ++delta_c_it) { *damage_it = std::max((*delta_max_it - delta_0) / (*delta_c_it - delta_0), 0.); *damage_it = std::min(*damage_it, 1.); } } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveBilinear); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc index c755e58f9..74a86753b 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc @@ -1,363 +1,362 @@ /** * @file material_cohesive_exponential.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Mon Jul 09 2012 * @date last modification: Fri Mar 21 2014 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @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 "material_cohesive_exponential.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveExponential::MaterialCohesiveExponential(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id) { AKANTU_DEBUG_IN(); this->registerParam("beta" , beta , 0. , _pat_parsable, "Beta parameter" ); - this->registerParam("delta_c", delta_c, 0. , _pat_parsable, "Critical displacement" ); this->registerParam("exponential_penalty", exp_penalty, true, _pat_parsable, "Is contact penalty following the exponential law?" ); this->registerParam("contact_tangent", contact_tangent, 1.0, _pat_parsable, "Ratio of contact tangent over the initial exponential tangent" ); // this->initInternalArray(delta_max, 1, _ek_cohesive); use_previous_delta_max=true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); if ((exp_penalty)&&(contact_tangent!=1)){ contact_tangent = 1; AKANTU_DEBUG_WARNING("The parsed paramter is forced to 1.0 when the contact penalty follows the exponential law"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::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); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator traction_end = tractions(el_type, ghost_type).end(spatial_dimension); Array::iterator delta_max_it = delta_max(el_type, ghost_type).begin(); Array::iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); /// compute scalars Real beta2 = beta*beta; /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++delta_max_it, ++delta_max_prev_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \beta^2 \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); if ((normal_opening_norm<0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { Vector op_n(*opening_it); op_n *= *normal_it; Vector delta_s(*opening_it); delta_s -= op_n; delta = tangential_opening_norm*beta; computeCoupledTraction(*traction_it, *normal_it, delta, delta_s, *delta_max_it, *delta_max_prev_it); computeCompressiveTraction(*traction_it, *normal_it, normal_opening_norm, *opening_it); } else computeCoupledTraction(*traction_it, *normal_it, delta, *opening_it, *delta_max_it, *delta_max_prev_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTraction(Vector & tract, const Vector & normal, Real delta, const Vector & opening, Real & delta_max_new, Real delta_max) { AKANTU_DEBUG_IN(); /// full damage case if (std::abs(delta) < Math::getTolerance()) { /// set traction to zero tract.clear(); } else { /// element not fully damaged /** * Compute traction loading @f$ \mathbf{T} = * e \sigma_c \frac{\delta}{\delta_c} e^{-\delta/ \delta_c}@f$ */ /** * Compute traction unloading @f$ \mathbf{T} = * \frac{t_{max}}{\delta_{max}} \delta @f$ */ Real beta2 = beta*beta; Vector op_n_n(opening); op_n_n*=normal; op_n_n*=(1-beta2); op_n_n*=normal; tract = beta2*opening; tract += op_n_n; delta_max_new = std::max(delta_max, delta); tract *= exp(1)*sigma_c*exp(- delta_max_new / delta_c)/delta_c; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressiveTraction(Vector & tract, const Vector & normal, Real delta_n, const Vector & opening) { Vector temp_tract(normal); if(exp_penalty) { temp_tract *= delta_n*exp(1)*sigma_c*exp(- delta_n / delta_c)/delta_c; } else { Real initial_tg = contact_tangent*exp(1)*sigma_c*delta_n/delta_c; temp_tract *= initial_tg; } tract += temp_tract; } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); 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); Array::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::iteratordelta_max_it = delta_max.previous(el_type, ghost_type).begin(); Real beta2 = beta*beta; /** * compute tangent matrix @f$ \frac{\partial \mathbf{t}} * {\partial \delta} = \hat{\mathbf{t}} \otimes * \frac{\partial (t/\delta)}{\partial \delta} * \frac{\hat{\mathbf{t}}}{\delta}+ \frac{t}{\delta} [ \beta^2 \mathbf{I} + * (1-\beta^2) (\mathbf{n} \otimes \mathbf{n})] @f$ **/ /** * In which @f$ * \frac{\partial(t/ \delta)}{\partial \delta} = * \left\{\begin{array} {l l} * -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad if * \delta \geq \delta_{max} \\ * 0 & \quad if \delta < \delta_{max}, \delta_n > 0 * \end{array}\right. @f$ **/ for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++traction_it, ++ delta_max_it) { Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); if ((normal_opening_norm<0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { Vector op_n(*opening_it); op_n *= *normal_it; Vector delta_s(*opening_it); delta_s -= op_n; delta = tangential_opening_norm*beta; computeCoupledTangent(*tangent_it, *traction_it, *normal_it, delta, delta_s, *delta_max_it); computeCompressivePenalty(*tangent_it, *normal_it, normal_opening_norm); } else computeCoupledTangent(*tangent_it, *traction_it, *normal_it, delta, *opening_it, *delta_max_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTangent(Matrix & tangent, Vector tract, const Vector & normal, Real delta, const Vector & opening, Real delta_max_new) { AKANTU_DEBUG_IN(); Real beta2 = beta*beta; Matrix I(spatial_dimension, spatial_dimension); I.eye(beta2); Matrix nn(spatial_dimension, spatial_dimension); nn.outerProduct(normal, normal); nn *= (1-beta2); I += nn; if(std::abs(delta) < Math::getTolerance()){ I *= exp(1)*sigma_c/delta_c; tangent += I; } else { Real traction_norm = tract.norm(); Vector t_hat(opening); Vector beta2delta(opening); beta2delta *= beta2; t_hat *= normal; t_hat *= (beta2-1); t_hat *= normal; t_hat += beta2delta; Vector deriv_t_hat(t_hat); deriv_t_hat /= delta; Real derivatives = 0; if ((delta > delta_max_new) || (std::abs(delta - delta_max_new) < 1e-12)) { derivatives = -exp(1- delta/delta_c)*sigma_c/(delta_c*delta_c); } deriv_t_hat *= derivatives; Matrix t_var_t(spatial_dimension, spatial_dimension); t_var_t.outerProduct(t_hat, deriv_t_hat); I *= traction_norm / delta; tangent += t_var_t; tangent += I; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressivePenalty(Matrix & tangent, const Vector & normal, Real delta_n) { if (!exp_penalty) delta_n=0; Matrix n_outer_n(spatial_dimension,spatial_dimension); n_outer_n.outerProduct(normal,normal); Real normal_tg = contact_tangent*exp(1)*sigma_c*exp(-delta_n/delta_c)*(1-delta_n/delta_c)/delta_c; n_outer_n *= normal_tg; tangent += n_outer_n; } INSTANSIATE_MATERIAL(MaterialCohesiveExponential); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh index 03a8d2fc0..6bfdb33d3 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh @@ -1,129 +1,126 @@ /** * @file material_cohesive_exponential.hh * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Mon Jul 09 2012 * @date last modification: Mon May 13 2013 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @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 "material_cohesive.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ #define __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material Exponential damage * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: 0) * - delta_c : critical opening (default: 0) */ template class MaterialCohesiveExponential : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveExponential(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// Initialization void initMaterial(); /// constitutive law void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type = _not_ghost); private: void computeCoupledTraction(Vector & tract, const Vector & normal, Real delta, const Vector & opening, Real & delta_max_new, Real delta_max); void computeCompressiveTraction(Vector & tract, const Vector & normal, Real delta_n, const Vector & opening); void computeCoupledTangent(Matrix & tangent, Vector tract, const Vector & normal, Real delta, const Vector & opening, Real delta_max_new); void computeCompressivePenalty(Matrix & tangent, const Vector & normal, Real delta_n); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// beta parameter Real beta; - /// critical displacement - Real delta_c; - /// contact penalty = initial slope ? bool exp_penalty; /// Ratio of contact tangent over the initial exponential tangent Real contact_tangent; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ // #include "material_cohesive_exponential_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ */ 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 78b4687ac..7d5c3ec56 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,503 +1,517 @@ /** * @file material_cohesive_linear.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Thu Aug 07 2014 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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 "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("delta_c", *this), + delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta" , beta , 0. , _pat_parsable | _pat_readable, "Beta parameter" ); - this->registerParam("G_cI" , G_cI , 0. , + + this->registerParam("G_c" , G_c , 0. , _pat_parsable | _pat_readable, "Mode I fracture energy" ); - this->registerParam("G_cII" , G_cII , 0. , - _pat_parsable | _pat_readable, - "Mode II fracture energy"); + this->registerParam("penalty", penalty, 0. , _pat_parsable | _pat_readable, "Penalty coefficient" ); + this->registerParam("volume_s", volume_s, 0. , _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); + this->registerParam("m_s", m_s, 1. , _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); - this->registerParam("kappa" , kappa , 1. , _pat_readable, "Kappa parameter"); + + this->registerParam("kappa" , kappa , 1. , + _pat_parsable | _pat_readable, + "Kappa parameter"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); - if (G_cII != 0) - kappa = G_cII / G_cI; - /// compute scalars beta2_kappa2 = beta * beta/kappa/kappa; beta2_kappa = beta * beta/kappa; - if (beta == 0) + if (Math::are_float_equal(beta, 0)) beta2_inv = 0; else beta2_inv = 1./beta/beta; - sigma_c_eff .initialize( 1); - delta_c .initialize( 1); + sigma_c_eff.initialize(1); + delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); + 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(); 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.getNbQuadraturePoints(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.getNbQuadraturePoints(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() { 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(type_cohesive); + Array & del_c = delta_c_eff(type_cohesive); Array & ins_stress = insertion_stress(type_cohesive); Array & trac_old = tractions_old(type_cohesive); const Array & f_stress = model->getStressOnFacets(type_facet); const Array & sigma_lim = sigma_c(type_facet); UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet); UInt nb_facet = f_filter.getSize(); if (nb_facet == 0) continue; UInt tot_nb_quad = nb_facet * nb_quad_facet; Array stress_check (tot_nb_quad); Array normal_traction(tot_nb_quad, spatial_dimension); computeStressNorms(f_stress, stress_check, normal_traction, type_facet); Array::iterator > stress_check_it = stress_check.begin_reinterpret(nb_quad_facet, nb_facet); 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 f = 0; f < nb_facet; ++f, ++sigma_lim_it, ++stress_check_it, ++normal_traction_it) { UInt facet = f_filter(f); if (!facets_check(facet)) continue; Real mean_stress = std::accumulate(stress_check_it->storage(), stress_check_it->storage() + nb_quad_facet, 0.); mean_stress /= nb_quad_facet; if (mean_stress > (*sigma_lim_it - tolerance)) { f_insertion(facet) = true; for (UInt q = 0; q < nb_quad_facet; ++q) { Real new_sigma = (*stress_check_it)(q); - Real new_delta = 2 * G_cI / new_sigma; Vector ins_s(normal_traction_it->storage() + q * spatial_dimension, spatial_dimension); if (spatial_dimension != 3) ins_s *= -1.; sig_c_eff.push_back(new_sigma); - del_c.push_back(new_delta); ins_stress.push_back(ins_s); trac_old.push_back(ins_s); + + 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; + + del_c.push_back(new_delta); } #if defined (AKANTU_DEBUG_TOOLS) && defined(AKANTU_CORE_CXX11) debug::element_manager.print(debug::_dm_material_cohesive, [facet, type_facet, nb_quad_facet, &f_stress](const Element & el)->std::string { std::stringstream sout; Element facet_el(type_facet, facet, _not_ghost); if(facet_el == el) { Array::const_iterator< Vector > stress = f_stress.begin(f_stress.getNbComponent()); stress += nb_quad_facet * facet; for (UInt qs = 0; qs < nb_quad_facet; ++qs, ++stress) { sout << *stress; } } return sout.str(); }); #endif } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_traction, Real & effective_norm) { AKANTU_DEBUG_IN(); normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (spatial_dimension == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (spatial_dimension == 3) { for (UInt s = 0; s < spatial_dimension - 1; ++s) { const Vector tangent_v(tangent.storage() + s * spatial_dimension, spatial_dimension); Real tangent_contrib_tmp = normal_traction.dot(tangent_v); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } } tangent_contrib = std::sqrt(tangent_contrib); normal_contrib = std::max(0., normal_contrib); effective_norm = std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeStressNorms(const Array & facet_stress, Array & stress_check, Array & normal_traction, ElementType type_facet) { AKANTU_DEBUG_IN(); Array & facets_check = model->getElementInserter().getCheckFacets(type_facet); Array & f_filter = facet_filter(type_facet); UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet); const Array & tangents = model->getTangents(type_facet); const Array & normals = model->getFEEngine("FacetsFEEngine").getNormalsOnQuadPoints(type_facet); Real * stress_check_it = stress_check.storage(); Array::const_iterator< Vector > normal_begin = normals.begin(spatial_dimension); Array::const_iterator< Vector > tangent_begin = tangents.begin(tangents.getNbComponent()); Array::const_iterator< Matrix > facet_stress_begin = facet_stress.begin(spatial_dimension, spatial_dimension * 2); Array::iterator > normal_traction_it = normal_traction.begin(spatial_dimension); Matrix stress_tmp(spatial_dimension, spatial_dimension); UInt nb_facet = f_filter.getSize(); UInt sp2 = spatial_dimension * spatial_dimension; UInt * current_facet = f_filter.storage(); for (UInt f = 0; f < nb_facet; ++f, ++current_facet) { if (!facets_check(*current_facet)) { stress_check_it += nb_quad_facet; normal_traction_it += nb_quad_facet; continue; } UInt current_quad = *current_facet * nb_quad_facet; for (UInt q = 0; q < nb_quad_facet; ++q, ++stress_check_it, ++normal_traction_it, ++current_quad) { 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 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.; /// compute normal, tangential and effective stress computeEffectiveNorm(stress_tmp, normal, tangent, *normal_traction_it, *stress_check_it); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::iterator< Vector > traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::const_iterator< Vector > normal_it = normal.begin(spatial_dimension); Array::iterator< Vector >traction_end = tractions(el_type, ghost_type).end(spatial_dimension); Array::iteratorsigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Array::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); Array::iteratordelta_c_it = - delta_c(el_type, ghost_type).begin(); + delta_c_eff(el_type, ghost_type).begin(); Array::iteratordamage_it = damage(el_type, ghost_type).begin(); Array::iterator > insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Real * memory_space = new Real[2*spatial_dimension]; Vector normal_opening(memory_space, spatial_dimension); Vector tangential_opening(memory_space + spatial_dimension, spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm * tangential_opening_norm * beta2_kappa2; bool penetration = normal_opening_norm < -Math::getTolerance(); if (penetration) { /// use penalty coefficient in case of penetration *contact_traction_it = normal_opening; *contact_traction_it *= penalty; *contact_opening_it = normal_opening; /// don't consider penetration contribution for delta *opening_it = tangential_opening; normal_opening.clear(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction_it->clear(); contact_opening_it->clear(); } delta = std::sqrt(delta); /// update maximum displacement and damage *delta_max_it = std::max(*delta_max_it, delta); *damage_it = std::min(*delta_max_it / *delta_c_it, 1.); /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ if (Math::are_float_equal(*damage_it, 1.)) traction_it->clear(); else if (Math::are_float_equal(*damage_it, 0.)) { if (penetration) traction_it->clear(); else *traction_it = *insertion_stress_it; } else { *traction_it = tangential_opening; *traction_it *= beta2_kappa; *traction_it += normal_opening; AKANTU_DEBUG_ASSERT(*delta_max_it != 0., "Division by zero, tolerance might be too low"); *traction_it *= *sigma_c_it / *delta_max_it * (1. - *damage_it); } } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveLinear); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh index 867fd039c..037fca00f 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh @@ -1,162 +1,160 @@ /** * @file material_cohesive_linear.hh * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Tue Jul 29 2014 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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 "material_cohesive.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material linear damage for extrinsic case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: 0) * - G_cI : fracture energy for mode I (default: 0) * - G_cII : fracture energy for mode II (default: 0) * - penalty : stiffness in compression to prevent penetration */ template class MaterialCohesiveLinear : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveLinear(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); /// check stress for cohesive elements' insertion virtual void checkInsertion(); protected: /// constitutive law void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /// compute stress norms on quadrature points for each facet for stress check virtual void computeStressNorms(const Array & facet_stress, Array & stress_check, Array & normal_stress, ElementType type_facet); /// compute effective stress norm for insertion check inline void computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_stress, Real & effective_norm); /** * Scale insertion traction sigma_c according to the volume of the * two elements surrounding a facet */ void scaleInsertionTraction(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get sigma_c_eff AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(InsertionTraction, sigma_c_eff, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// beta parameter Real beta; /// beta square inverse to compute effective norm Real beta2_inv; /// mode I fracture energy - Real G_cI; - - /// mode II fracture energy - Real G_cII; + Real G_c; /// kappa parameter Real kappa; /// constitutive law scalar to compute delta Real beta2_kappa2; /// constitutive law scalar to compute traction Real beta2_kappa; /// penalty coefficient Real penalty; /// reference volume used to scale sigma_c Real volume_s; /// weibull exponent used to scale sigma_c Real m_s; /// critical effective stress RandomInternalField sigma_c_eff; - /// critical displacement - CohesiveInternalField delta_c; + /// effective critical displacement (each element can have a + /// different value) + CohesiveInternalField delta_c_eff; /// stress at insertion CohesiveInternalField insertion_stress; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_cohesive_linear_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc index 14376a897..184d43a1f 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,616 +1,618 @@ /** * @file material_cohesive.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jul 29 2014 * * @brief Specialization of the material class 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 "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" #include "aka_random_generator.hh" #include "shape_cohesive.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive(&(model.getFEEngineClass("CohesiveFEEngine"))), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), opening_old("opening (old)", *this), tractions("tractions", *this), tractions_old("tractions (old)", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta max", *this), use_previous_delta_max(false), damage("damage", *this), sigma_c("sigma_c", *this) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); + this->registerParam("delta_c", delta_c, 0., + _pat_parsable, "Critical displacement" ); this->model->getMesh().initElementTypeMapArray(this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); this->reversible_energy.initialize(1 ); this->total_energy .initialize(1 ); this->tractions_old .initialize(spatial_dimension); this->tractions .initialize(spatial_dimension); this->opening_old .initialize(spatial_dimension); this->contact_tractions.initialize(spatial_dimension); this->contact_opening .initialize(spatial_dimension); this->opening .initialize(spatial_dimension); this->delta_max .initialize(1 ); this->damage .initialize(1 ); if (model->getIsExtrinsic()) this->sigma_c.initialize(1); if (use_previous_delta_max) delta_max.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif Array & residual = const_cast &>(model->getResidual()); 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); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & traction = tractions(*it, ghost_type); Array & contact_traction = contact_tractions(*it, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); /// compute @f$t_i N_a@f$ Array * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); Array::iterator > traction_it = traction.begin(spatial_dimension, 1); Array::iterator > contact_traction_it = contact_traction.begin(spatial_dimension, 1); Array::const_iterator > shapes_filtered_begin = shapes.begin(1, size_of_shapes); Array::iterator > traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ Array * int_t_N = new Array(nb_element, spatial_dimension*size_of_shapes, "int_t_N"); fem_cohesive->integrate(*traction_cpy, *int_t_N, spatial_dimension * size_of_shapes, *it, ghost_type, elem_filter); delete traction_cpy; int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent()); Real * int_t_N_val = int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element * spatial_dimension; } /// assemble model->getFEEngineBoundary().assembleArray(*int_t_N, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, 1); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); 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); for(; it != last_type; ++it) { UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if(!nb_element) continue; UInt size_of_shapes = shapes.getNbComponent(); Array * shapes_filtered = new Array(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_val = shapes.storage(); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } /** * compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c c c} * 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ * 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ * 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ * 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ * 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\ * 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 * \end{array} \right]@f$ **/ // UInt size_of_A = spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element; // Real * A = new Real[size_of_A]; // memset(A, 0, size_of_A*sizeof(Real)); Matrix A(spatial_dimension*size_of_shapes, spatial_dimension*nb_nodes_per_element); for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension*size_of_shapes) = -1; } /// compute traction computeTraction(ghost_type); /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} @f$ Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * nb_quadrature_points, spatial_dimension, "normal"); Array normal(nb_quadrature_points, spatial_dimension, "normal"); computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); tangent_stiffness_matrix->clear(); computeTangentTraction(*it, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension*nb_nodes_per_element*spatial_dimension*nb_nodes_per_element; Array * at_nt_d_n_a = new Array (nb_element*nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator > shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N (spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A (spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for(; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.clear(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension ; ++i) for (UInt n = 0; n < size_of_shapes; ++n) N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} {\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; Array * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive->integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, *it, ghost_type, elem_filter); delete at_nt_d_n_a; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, *it, ghost_type, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif 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); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = nb_element*fem_cohesive->getNbQuadraturePoints(*it, ghost_type); Array normal(nb_quadrature_points, spatial_dimension, "normal"); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (type == _cohesive_1d_2) fem_cohesive->computeNormalsOnControlPoints(position, normal, type, ghost_type); else { #define COMPUTE_NORMAL(type) \ fem_cohesive->getShapeFunctions(). \ computeNormalsOnControlPoints(position, \ normal, \ ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_OPENING(type) \ fem_cohesive->getShapeFunctions(). \ interpolateOnControlPoints(displacement, \ opening, \ spatial_dimension, \ ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); Real * memory_space = new Real[2*spatial_dimension]; Vector b(memory_space, spatial_dimension); Vector h(memory_space + spatial_dimension, spatial_dimension); for(; it != last_type; ++it) { Array::iterator erev = reversible_energy(*it, _not_ghost).begin(); Array::iterator etot = total_energy(*it, _not_ghost).begin(); Array::vector_iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_end = tractions(*it, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } delete [] memory_space; /// update old values it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); GhostType ghost_type = _not_ghost; for(; it != last_type; ++it) { tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type)); opening_old(*it, ghost_type).copy(opening(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Array dissipated_energy(total_energy(*it, _not_ghost)); dissipated_energy -= reversible_energy(*it, _not_ghost); edis += fem_cohesive->integrate(dissipated_energy, *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Array & el_filter = element_filter(*it, _not_ghost); UInt nb_quad_per_el = fem_cohesive->getNbQuadraturePoints(*it, _not_ghost); UInt nb_quad_points = el_filter.getSize() * nb_quad_per_el; Array contact_energy(nb_quad_points); Array::vector_iterator contact_traction_it = contact_tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(*it, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt el = 0; el < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++el) { contact_energy(el) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive->integrate(contact_energy, *it, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); else if (type == "cohesive contact") return getContactEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh index 38600f1fb..2294ebec7 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh @@ -1,243 +1,246 @@ /** * @file material_cohesive.hh * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jul 29 2014 * * @brief Specialization of the material class 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 "aka_common.hh" #include "material.hh" #include "fe_engine_template.hh" #include "aka_common.hh" #include "cohesive_internal_field.hh" #include "cohesive_element_inserter.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_HH__ #define __AKANTU_MATERIAL_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } __BEGIN_AKANTU__ class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate MyFEEngineCohesiveType; public: MaterialCohesive(SolidMechanicsModel& model, const ID & id = ""); virtual ~MaterialCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleResidual(GhostType ghost_type = _not_ghost); /// compute reversible and total energies by element void computeEnergies(); /// check stress for cohesive elements' insertion virtual void checkInsertion() { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(__attribute__((unused)) const ElementType type, __attribute__((unused)) Array & result) { }; virtual void computeAllStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) { }; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) const Array & normal, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } void computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type); void computeOpening(const Array & displacement, Array & normal, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type); /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy virtual Real getEnergy(std::string type); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model MyFEEngineCohesiveType * fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// opening in all elements and quadrature points (previous time step) CohesiveInternalField opening_old; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction in all elements and quadrature points (previous time step) CohesiveInternalField tractions_old; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; + + /// critical displacement + Real delta_c; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" __END_AKANTU__ #include "cohesive_internal_field_tmpl.hh" #endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_1d_element/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_1d_element/material.dat index 42c6fb6fa..5f1e159c1 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_1d_element/material.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_1d_element/material.dat @@ -1,11 +1,11 @@ material elastic [ name = steel rho = 1e3 # density E = 1.e3 # young's modulus ] material cohesive_linear [ name = cohesive - G_cI = 100 + G_c = 100 sigma_c = 100 ] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/material.dat index 2c816e365..bfff2e8c3 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/material.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/material.dat @@ -1,16 +1,15 @@ seed = 1 material elastic [ name = bulk rho = 2500 E = 42.e9 nu = 0.2 ] material cohesive_linear [ name = cohesive beta = 1 - G_cI = 100 - G_cII = 100 + G_c = 100 sigma_c = 300e6 uniform [0, 30e6] ] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat index dcf925bab..7fceb2d23 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat @@ -1,23 +1,21 @@ material elastic [ name = steel rho = 1e3 # density E = 1e3 # young's modulus nu = 0.001 # poisson's ratio ] material cohesive_linear [ name = cohesive beta = 1 - G_cI = 100 - G_cII = 100 + G_c = 100 sigma_c = 100 volume_s = 1 ] material cohesive_linear [ name = cohesive2 beta = 1 - G_cI = 100 - G_cII = 100 + G_c = 100 sigma_c = 1000 ] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat index 0f9f4b536..d31ad403d 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat @@ -1,16 +1,15 @@ material elastic [ name = steel rho = 10 # density E = 10 # young's modulus nu = 0.1 # poisson's ratio ] material cohesive_bilinear [ name = cohesive sigma_c = 1 uniform [0, 1] beta = 1.5 - G_cI = 1 - G_cII = 1 + G_c = 1 delta_0 = 1 # penalty = 1e10 ] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material_tetrahedron.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material_tetrahedron.dat index 11213337d..ba242d682 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material_tetrahedron.dat +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material_tetrahedron.dat @@ -1,16 +1,16 @@ material elastic [ name = steel rho = 10 # density E = 10 # young's modulus nu = 0.1 # poisson's ratio ] material cohesive_bilinear [ name = cohesive sigma_c = 1 beta = 1.5 - G_cI = 1 - G_cII = 1.2 + G_c = 1 + kappa = 1.2 delta_0 = 0.1 # penalty = 1e10 ] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc index 6e7ca68ce..30d64c941 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic_tetrahedron.cc @@ -1,441 +1,440 @@ /** * @file test_cohesive_intrinsic_tetrahedron.cc * * @author Marco Vocialta * * @date creation: Tue Aug 27 2013 * @date last modification: Fri Sep 19 2014 * * @brief Test 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 #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" #include "material_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; void updateDisplacement(SolidMechanicsModelCohesive &, Array &, ElementType, Vector &); bool checkTractions(SolidMechanicsModelCohesive & model, ElementType type, Vector & opening, Vector & theoretical_traction, Matrix & rotation); void findNodesToCheck(const Mesh & mesh, const Array & elements, ElementType type, Array & nodes_to_check); bool checkEquilibrium(const Array & residual); bool checkResidual(const Array & residual, const Vector & traction, const Array & nodes_to_check, const Matrix & rotation); int main(int argc, char *argv[]) { initialize("material_tetrahedron.dat", argc, argv); // debug::setDebugLevel(dblDump); const UInt spatial_dimension = 3; const UInt max_steps = 60; const Real increment_constant = 0.01; Math::setTolerance(1.e-12); const ElementType type = _tetrahedron_10; Mesh mesh(spatial_dimension); mesh.read("tetrahedron.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(); model.limitInsertion(_x, -0.01, 0.01); model.insertIntrinsicElements(); Array & boundary = model.getBlockedDOFs(); boundary.set(true); UInt nb_element = mesh.getNbElement(type); model.updateResidual(); model.setBaseName("intrinsic_tetrahedron"); model.addDumpFieldVector("displacement"); model.addDumpField("residual"); model.dump(); model.setBaseNameToDumper("cohesive elements", "cohesive_elements_tetrahedron"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.dump("cohesive elements"); /// find elements to displace Array elements; Real * bary = new Real[spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter(el, type, bary); if (bary[0] > 0.01) elements.push_back(el); } delete[] bary; /// find nodes to check Array nodes_to_check; findNodesToCheck(mesh, elements, type, nodes_to_check); /// rotate mesh Real angle = 1.; Matrix rotation(spatial_dimension, spatial_dimension); rotation.clear(); rotation(0, 0) = std::cos(angle); rotation(0, 1) = std::sin(angle) * -1.; rotation(1, 0) = std::sin(angle); rotation(1, 1) = std::cos(angle); rotation(2, 2) = 1.; Vector increment_tmp(spatial_dimension); for (UInt dim = 0; dim < spatial_dimension; ++dim) { increment_tmp(dim) = (dim + 1) * increment_constant; } Vector increment(spatial_dimension); increment.mul(rotation, increment_tmp); Array & position = mesh.getNodes(); Array position_tmp(position); Array::iterator > position_it = position.begin(spatial_dimension); Array::iterator > position_end = position.end(spatial_dimension); Array::iterator > position_tmp_it = position_tmp.begin(spatial_dimension); for (; position_it != position_end; ++position_it, ++position_tmp_it) position_it->mul(rotation, *position_tmp_it); model.dump(); model.dump("cohesive elements"); updateDisplacement(model, elements, type, increment); Real theoretical_Ed = 0; Vector opening(spatial_dimension); Vector traction(spatial_dimension); Vector opening_old(spatial_dimension); Vector traction_old(spatial_dimension); opening.clear(); traction.clear(); opening_old.clear(); traction_old.clear(); Vector Dt(spatial_dimension); Vector Do(spatial_dimension); const Array & residual = model.getResidual(); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.updateResidual(); opening += increment_tmp; if (checkTractions(model, type, opening, traction, rotation) || checkEquilibrium(residual) || checkResidual(residual, traction, nodes_to_check, rotation)) { finalize(); return EXIT_FAILURE; } /// compute energy Do = opening; Do -= opening_old; Dt = traction_old; Dt += traction; theoretical_Ed += .5 * Do.dot(Dt); opening_old = opening; traction_old = traction; updateDisplacement(model, elements, type, increment); if(s % 10 == 0) { std::cout << "passing step " << s << "/" << max_steps << std::endl; model.dump(); model.dump("cohesive elements"); } } model.dump(); model.dump("cohesive elements"); Real Ed = model.getEnergy("dissipated"); theoretical_Ed *= 4.; std::cout << Ed << " " << theoretical_Ed << std::endl; if (!Math::are_float_equal(Ed, theoretical_Ed) || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic_tetrahedron was passed!" << std::endl; return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void updateDisplacement(SolidMechanicsModelCohesive & model, Array & elements, ElementType type, Vector & increment) { UInt spatial_dimension = model.getSpatialDimension(); Mesh & mesh = model.getFEEngine().getMesh(); UInt nb_element = elements.getSize(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Array & connectivity = mesh.getConnectivity(type); Array & displacement = model.getDisplacement(); Array update(nb_nodes); update.clear(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity(elements(el), n); if (!update(node)) { Vector node_disp(displacement.storage() + node * spatial_dimension, spatial_dimension); node_disp += increment; update(node) = true; } } } } /* -------------------------------------------------------------------------- */ bool checkTractions(SolidMechanicsModelCohesive & model, ElementType type, Vector & opening, Vector & theoretical_traction, Matrix & rotation) { UInt spatial_dimension = model.getSpatialDimension(); const MaterialCohesive & mat_cohesive = dynamic_cast < const MaterialCohesive & > (model.getMaterial(1)); Real sigma_c = mat_cohesive.getParam< RandomInternalField >("sigma_c"); const Real beta = mat_cohesive.getParam("beta"); - const Real G_cI = mat_cohesive.getParam("G_cI"); - // Real G_cII = mat_cohesive.getParam("G_cII"); + const Real G_c = mat_cohesive.getParam("G_c"); const Real delta_0 = mat_cohesive.getParam("delta_0"); const Real kappa = mat_cohesive.getParam("kappa"); - Real delta_c = 2 * G_cI / sigma_c; + Real delta_c = 2 * G_c / sigma_c; sigma_c *= delta_c / (delta_c - delta_0); ElementType type_facet = Mesh::getFacetType(type); ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const Array & traction = mat_cohesive.getTraction(type_cohesive); const Array & damage = mat_cohesive.getDamage(type_cohesive); UInt nb_quad_per_el = model.getFEEngine("CohesiveFEEngine").getNbQuadraturePoints(type_cohesive); UInt nb_element = model.getMesh().getNbElement(type_cohesive); UInt tot_nb_quad = nb_element * nb_quad_per_el; Vector normal_opening(spatial_dimension); normal_opening.clear(); normal_opening(0) = opening(0); Real normal_opening_norm = normal_opening.norm(); Vector tangential_opening(spatial_dimension); tangential_opening.clear(); for (UInt dim = 1; dim < spatial_dimension; ++dim) tangential_opening(dim) = opening(dim); Real tangential_opening_norm = tangential_opening.norm(); Real beta2_kappa2 = beta * beta/kappa/kappa; Real beta2_kappa = beta * beta/kappa; Real delta = std::sqrt(tangential_opening_norm * tangential_opening_norm * beta2_kappa2 + normal_opening_norm * normal_opening_norm); delta = std::max(delta, delta_0); Real theoretical_damage = std::min(delta / delta_c, 1.); if (Math::are_float_equal(theoretical_damage, 1.)) theoretical_traction.clear(); else { theoretical_traction = tangential_opening; theoretical_traction *= beta2_kappa; theoretical_traction += normal_opening; theoretical_traction *= sigma_c / delta * (1. - theoretical_damage); } // adjust damage theoretical_damage = std::max((delta - delta_0) / (delta_c - delta_0), 0.); theoretical_damage = std::min(theoretical_damage, 1.); Vector theoretical_traction_rotated(spatial_dimension); theoretical_traction_rotated.mul(rotation, theoretical_traction); for (UInt q = 0; q < tot_nb_quad; ++q) { for (UInt dim = 0; dim < spatial_dimension; ++dim) { if (!Math::are_float_equal(theoretical_traction_rotated(dim), traction(q, dim))) { std::cout << "Tractions are incorrect" << std::endl; return 1; } } if (!Math::are_float_equal(theoretical_damage, damage(q))) { std::cout << "Damage is incorrect" << std::endl; return 1; } } return 0; } /* -------------------------------------------------------------------------- */ void findNodesToCheck(const Mesh & mesh, const Array & elements, ElementType type, Array & nodes_to_check) { const Array & connectivity = mesh.getConnectivity(type); const Array & position = mesh.getNodes(); UInt nb_nodes = position.getSize(); UInt nb_nodes_per_elem = connectivity.getNbComponent(); Array checked_nodes(nb_nodes); checked_nodes.clear(); for (UInt el = 0; el < elements.getSize(); ++el) { UInt element = elements(el); Vector conn_el(connectivity.storage() + nb_nodes_per_elem * element, nb_nodes_per_elem); for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = conn_el(n); if (Math::are_float_equal(position(node, 0), 0.) && checked_nodes(node) == false) { checked_nodes(node) = true; nodes_to_check.push_back(node); } } } } /* -------------------------------------------------------------------------- */ bool checkEquilibrium(const Array & residual) { UInt spatial_dimension = residual.getNbComponent(); Vector residual_sum(spatial_dimension); residual_sum.clear(); Array::const_iterator > res_it = residual.begin(spatial_dimension); Array::const_iterator > res_end = residual.end(spatial_dimension); for (; res_it != res_end; ++res_it) residual_sum += *res_it; for (UInt s = 0; s < spatial_dimension; ++s) { if (!Math::are_float_equal(residual_sum(s), 0.)) { std::cout << "System is not in equilibrium!" << std::endl; return 1; } } return 0; } /* -------------------------------------------------------------------------- */ bool checkResidual(const Array & residual, const Vector & traction, const Array & nodes_to_check, const Matrix & rotation) { UInt spatial_dimension = residual.getNbComponent(); Vector total_force(spatial_dimension); total_force.clear(); for (UInt n = 0; n < nodes_to_check.getSize(); ++n) { UInt node = nodes_to_check(n); Vector res(residual.storage() + node * spatial_dimension, spatial_dimension); total_force += res; } Vector theoretical_total_force(spatial_dimension); theoretical_total_force.mul(rotation, traction); theoretical_total_force *= -1 * 2 * 2; for (UInt s = 0; s < spatial_dimension; ++s) { if (!Math::are_float_equal(total_force(s), theoretical_total_force(s))) { std::cout << "Total force isn't correct!" << std::endl; return 1; } } return 0; }