diff --git a/packages/00_core.cmake b/packages/00_core.cmake index 851e5afc1..4d0ec8a0c 100644 --- a/packages/00_core.cmake +++ b/packages/00_core.cmake @@ -1,349 +1,351 @@ #=============================================================================== # @file core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date Mon Nov 21 18:19:15 2011 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== set(AKANTU_CORE ON CACHE INTERNAL "core package for Akantu" FORCE) set(AKANTU_CORE_FILES common/aka_blas_lapack.hh common/aka_circular_vector.hh common/aka_circular_vector_inline_impl.cc common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.cc common/aka_csr.hh common/aka_error.cc common/aka_error.hh common/aka_event_handler.hh common/aka_extern.cc common/aka_fwd.hh common/aka_grid.hh common/aka_grid_dynamic.hh common/aka_grid_tmpl.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_memory.cc common/aka_memory.hh common/aka_memory_inline_impl.cc common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_static_memory.cc common/aka_static_memory.hh common/aka_static_memory_inline_impl.cc common/aka_static_memory_tmpl.hh common/aka_typelist.hh common/aka_types.hh common/aka_vector.cc common/aka_vector.hh common/aka_vector_tmpl.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper.cc fem/by_element_type.hh fem/by_element_type_tmpl.hh fem/element_class.cc fem/element_class.hh fem/element_class_tmpl.hh fem/element_classes/element_class_hexahedron_8_inline_impl.cc fem/element_classes/element_class_pentahedron_6_inline_impl.cc fem/element_classes/element_class_point_1_inline_impl.cc fem/element_classes/element_class_quadrangle_4_inline_impl.cc fem/element_classes/element_class_quadrangle_8_inline_impl.cc fem/element_classes/element_class_segment_2_inline_impl.cc fem/element_classes/element_class_segment_3_inline_impl.cc fem/element_classes/element_class_tetrahedron_10_inline_impl.cc fem/element_classes/element_class_tetrahedron_4_inline_impl.cc fem/element_classes/element_class_triangle_3_inline_impl.cc fem/element_classes/element_class_triangle_6_inline_impl.cc fem/element_group.cc fem/element_group.hh fem/element_group_inline_impl.cc fem/fem.cc fem/fem.hh fem/fem_inline_impl.cc fem/fem_template.hh fem/fem_template_tmpl.hh fem/geometrical_data_tmpl.hh fem/geometrical_element.cc fem/group_manager.cc fem/group_manager.hh fem/group_manager_inline_impl.cc fem/integration_element.cc fem/integrator.hh fem/integrator_gauss.hh fem/integrator_gauss_inline_impl.cc fem/interpolation_element.cc fem/interpolation_element_tmpl.hh fem/mesh.cc fem/mesh.hh fem/mesh_filter.hh fem/mesh_data.cc fem/mesh_data.hh fem/mesh_data_tmpl.hh fem/mesh_inline_impl.cc fem/node_group.cc fem/node_group.hh fem/node_group_inline_impl.cc fem/shape_functions.hh fem/shape_functions_inline_impl.cc fem/shape_lagrange.cc fem/shape_lagrange.hh fem/shape_lagrange_inline_impl.cc fem/shape_linked.cc fem/shape_linked.hh fem/shape_linked_inline_impl.cc io/dumper/dumpable.hh io/dumper/dumpable_inline_impl.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh io/model_io.cc io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parsable_tmpl.hh io/parser/parser.cc io/parser/parser.hh io/parser/parser_tmpl.hh mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.cc model/boundary_condition.hh model/boundary_condition_functor.hh model/boundary_condition_functor_inline_impl.cc model/boundary_condition_tmpl.hh model/integration_scheme/generalized_trapezoidal.hh model/integration_scheme/generalized_trapezoidal_inline_impl.cc model/integration_scheme/integration_scheme_1st_order.hh model/integration_scheme/integration_scheme_2nd_order.hh model/integration_scheme/newmark-beta.hh model/integration_scheme/newmark-beta_inline_impl.cc model/model.cc model/model.hh model/model_inline_impl.cc model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.cc model/solid_mechanics/material_list.hh model/solid_mechanics/material_random_internal.hh model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh solver/solver.cc solver/solver.hh solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_inline_impl.cc synchronizer/communication_buffer.hh synchronizer/communication_buffer_inline_impl.cc synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/data_accessor_inline_impl.cc synchronizer/distributed_synchronizer.cc synchronizer/distributed_synchronizer.hh synchronizer/distributed_synchronizer_tmpl.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.cc synchronizer/filtered_synchronizer.cc synchronizer/filtered_synchronizer.hh synchronizer/mpi_type_wrapper.hh synchronizer/pbc_synchronizer.cc synchronizer/pbc_synchronizer.hh synchronizer/real_static_communicator.hh synchronizer/static_communicator.cc synchronizer/static_communicator.hh synchronizer/static_communicator_dummy.hh synchronizer/static_communicator_inline_impl.hh synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh ) set(AKANTU_CORE_DEB_DEPEND libboost-dev ) set(AKANTU_CORE_TESTS test_csr test_facet_element_mapping test_facet_extraction_tetrahedron_4 test_facet_extraction_triangle_3 test_grid test_interpolate_stress test_local_material test_material_damage_non_local test_material_thermal test_matrix test_mesh_boundary test_mesh_data test_mesh_io_msh test_mesh_io_msh_physical_names test_mesh_partitionate_mesh_data test_parser test_pbc_tweak test_purify_mesh test_solid_mechanics_model_bar_traction2d test_solid_mechanics_model_bar_traction2d_structured test_solid_mechanics_model_bar_traction2d_structured_pbc test_solid_mechanics_model_boundary_condition test_solid_mechanics_model_circle_2 test_solid_mechanics_model_cube3d test_solid_mechanics_model_cube3d_pbc test_solid_mechanics_model_cube3d_tetra10 test_solid_mechanics_model_square + test_solid_mechanics_model_reassign_material + test_solid_mechanics_model_prestrain_material test_static_memory test_surface_extraction_tetrahedron_4 test_surface_extraction_triangle_3 test_vector test_vector_iterator test_weight ) set(AKANTU_CORE_MANUAL_FILES manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-solidmechanicsmodel.tex manual-lumping.tex manual-elements.tex manual-appendix-elements.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg figures/doc_wheel.pdf figures/doc_wheel.svg figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/hooke_law.pdf figures/hot-point-1.png figures/hot-point-2.png figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/implicit_static.pdf figures/implicit_static.svg figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/law.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/vectors.pdf figures/vectors.svg figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) include(CheckFunctionExists) check_function_exists(clock_gettime _clock_gettime) if(NOT _clock_gettime) set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON) else() set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF) endif() diff --git a/packages/30_extra_materials.cmake b/packages/30_extra_materials.cmake index 88308f2f4..21d3585b6 100644 --- a/packages/30_extra_materials.cmake +++ b/packages/30_extra_materials.cmake @@ -1,98 +1,99 @@ #=============================================================================== # @file extra_materials.cmake # # @author Nicolas Richart # # @date Wed Oct 31 16:24:42 2012 # # @brief package description for extra materials list # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== option(AKANTU_EXTRA_MATERIALS "Add the extra list of materials in Akantu" OFF) add_external_package_dependencies(extra_materials lapack) set(AKANTU_EXTRA_MATERIALS_FILES model/solid_mechanics/materials/material_neohookean.cc model/solid_mechanics/materials/material_viscoelastic/material_stiffness_proportional.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_damage_linear.cc model/solid_mechanics/materials/material_extra_includes.hh model/solid_mechanics/materials/material_viscoelastic/material_stiffness_proportional.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_neohookean.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_damage_linear.hh model/solid_mechanics/materials/material_damage/material_vreepeerlings.hh model/solid_mechanics/materials/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_damage/material_damage_linear_inline_impl.cc model/solid_mechanics/materials/material_damage/material_vreepeerlings_inline_impl.cc model/solid_mechanics/materials/material_damage/material_vreepeerlings_tmpl.hh model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/solid_mechanics/materials/material_finite_deformation.cc model/solid_mechanics/materials/material_finite_deformation/material_plastic.cc model/solid_mechanics/materials/material_finite_deformation.hh model/solid_mechanics/materials/material_finite_deformation/material_plastic.hh model/solid_mechanics/materials/material_finite_deformation/material_plastic_inline_impl.cc model/solid_mechanics/materials/material_finite_deformation_inline_impl.cc model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc.hh model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc.cc model/solid_mechanics/materials/material_plasticityinc/material_plasticityinc_inline_impl.cc model/solid_mechanics/materials/material_viscoplasticity/material_viscoplasticity.hh model/solid_mechanics/materials/material_viscoplasticity/material_viscoplasticity.cc model/solid_mechanics/materials/material_viscoplasticity/material_viscoplasticity_inline_impl.cc ) set(AKANTU_EXTRA_MATERIALS_TESTS test_material_standard_linear_solid_deviatoric_relaxation test_material_standard_linear_solid_deviatoric_relaxation_tension test_material_plasticity + test_material_damage ) set(AKANTU_EXTRA_MATERIALS_DOC manual/manual-extra_materials.tex ) set(AKANTU_EXTRA_MATERIALS_MANUAL_FILES manual-extra_materials.tex figures/stress_strain_neo.pdf figures/stress_strain_visco.pdf figures/visco_elastic_law.pdf ) \ No newline at end of file diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 4a6575e73..ce03a0bce 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1358 +1,1526 @@ /** * @file material.cc * * @author Nicolas Richart * @author Marco Vocialta * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), finite_deformation(false), inelastic_deformation(false), name(""), model(&model), spatial_dimension(this->model->getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), + pre_strain("pre_strain", *this), strain("strain", *this), delta_stress("delta_stress", *this), delta_strain("delta_strain", *this), inelas_strain("inelas_strain", *this), piola_kirchhoff_stress("piola_kirchhoff_stress", *this), // potential_energy_vector(false), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), previous_stress("previous_stress", *this), use_previous_strain(false), previous_strain("previous_strain", *this), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the material model.getMesh().initByElementTypeArray(element_filter, 1, spatial_dimension, false, _ek_regular); registerParam("rho" , rho , 0. , _pat_parsable | _pat_modifiable, "Density"); registerParam("name" , name , std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation" , finite_deformation , false , _pat_parsable | _pat_modifiable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false , _pat_parsable | _pat_modifiable, "Is inelastic deformation"); /// allocate strain stress for local elements + pre_strain.initialize(spatial_dimension * spatial_dimension); strain.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if(finite_deformation) { this->delta_strain.initialize(spatial_dimension * spatial_dimension); this->delta_stress.initialize(spatial_dimension * spatial_dimension); this->piola_kirchhoff_stress.initialize(spatial_dimension * spatial_dimension); } if(inelastic_deformation) { this->delta_strain.initialize(spatial_dimension * spatial_dimension); this->delta_stress.initialize(spatial_dimension * spatial_dimension); this->inelas_strain.initialize(spatial_dimension * spatial_dimension); } if(use_previous_stress) { this->previous_stress.initialize(spatial_dimension * spatial_dimension); } if(use_previous_strain) { this->previous_strain.initialize(spatial_dimension * spatial_dimension); } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState(const GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getMesh().lastType(spatial_dimension, ghost_type); if(finite_deformation) for (; it != last_type; ++it) { if(use_previous_stress) previous_stress(*it, ghost_type).copy(piola_kirchhoff_stress(*it, ghost_type)); if(use_previous_strain) previous_strain(*it, ghost_type).copy(strain(*it, ghost_type)); } else for (; it != last_type; ++it) { if(use_previous_stress) previous_stress(*it, ghost_type).copy(stress(*it, ghost_type)); if(use_previous_strain) previous_strain(*it, ghost_type).copy(strain(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeAllStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); if(!finite_deformation){ Array & residual = const_cast &>(model->getResidual()); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Array & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Array * sigma_dphi_dx = new Array(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEM::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array & stress_vect = stress(*it, ghost_type); Array::iterator< Matrix > sigma = stress_vect.begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > B = shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::iterator< Matrix > Bt_sigma_it = sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element); for (UInt q = 0; q < nb_element*nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it) Bt_sigma_it->mul(*sigma, *B); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Array * int_sigma_dphi_dx = new Array(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model->getFEM().assembleArray(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete int_sigma_dphi_dx; } } else{ switch (spatial_dimension){ case 1: this->assembleResidual<1>(ghost_type); break; case 2: this->assembleResidual<2>(ghost_type); break; case 3: this->assembleResidual<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the strain * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); Array & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect, spatial_dimension, *it, ghost_type, elem_filter); + strain_vect += pre_strain(*it, ghost_type); + if(finite_deformation || inelastic_deformation){ /// compute @f$\nabla \delta u@f$ Array & delta_strain_vect = delta_strain(*it, ghost_type); model->getFEM().gradientOnQuadraturePoints(model->getIncrement(), delta_strain_vect, spatial_dimension, *it, ghost_type, elem_filter); } /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation,"The Cauchy stress can only be computed if you are working in finite deformation."); //resizeInternalArray(stress); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) switch (spatial_dimension){ case 1: this->computeCauchyStress<1>(*it, ghost_type); break; case 2: this->computeCauchyStress<2>(*it, ghost_type); break; case 3: this->computeCauchyStress<3>(*it, ghost_type); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::iterator< Matrix > strain_it = this->strain(el_type, ghost_type).begin(dim, dim); Array::iterator< Matrix > strain_end = this->strain(el_type, ghost_type).end(dim, dim); Array::iterator< Matrix > piola_it = this->piola_kirchhoff_stress(el_type, ghost_type).begin(dim, dim); Array::iterator< Matrix > stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; strain_it != strain_end; ++strain_it, ++piola_it, ++stress_it) { Matrix & grad_u = *strain_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF (grad_u, F_tensor); computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model->getDisplacement(); //resizeInternalArray(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); Array & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect, spatial_dimension, *it, ghost_type, elem_filter); setToSteadyState(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { if(finite_deformation){ switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL < 1 > (*it, ghost_type); assembleStiffnessMatrixL2 < 1 > (*it, ghost_type); break; } case 2: { assembleStiffnessMatrixNL < 2 > (*it, ghost_type); assembleStiffnessMatrixL2 < 2 > (*it, ghost_type); break; } case 3: { assembleStiffnessMatrixNL < 3 > (*it, ghost_type); assembleStiffnessMatrixL2 < 3 > (*it, ghost_type); break; } } } else { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(*it, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Array & shapes_derivatives = model->getFEM().getShapesDerivatives( type,ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & strain_vect = strain(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect, dim, type, ghost_type, elem_filter); if(inelastic_deformation){ /// compute @f$\nabla \delta u@f$ Array & delta_strain_vect = delta_strain(type, ghost_type); delta_strain_vect.resize(nb_quadrature_points * nb_element); model->getFEM().gradientOnQuadraturePoints(model->getIncrement(), delta_strain_vect, dim, type, ghost_type, elem_filter); } UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapesd_filtered = new Array(0, dim * nb_nodes_per_element, "filtered shapesd"); FEM::filterElementalData(model->getFEM().getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); Array::iterator< Matrix > shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::iterator< Matrix > Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim*nb_nodes_per_element); Array::iterator< Matrix > D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::iterator< Matrix > D_end = tangent_stiffness_matrix->end (tangent_size, tangent_size); for(; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) { Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEM().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); const Array & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); //Array & strain_vect = delta_strain(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); //strain_vect.resize(nb_quadrature_points * nb_element); // model->getFEM().gradientOnQuadraturePoints(model->getIncrement(), strain_vect, // dim, type, ghost_type, &elem_filter); Array * shapes_derivatives_filtered = new Array (nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_iterator< Matrix > shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::iterator< Matrix > shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array (nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::iterator< Matrix > Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); Array::iterator< Matrix > Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); Array::iterator< Matrix > piola_it = piola_kirchhoff_stress(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { Matrix & Bt_S_B = *Bt_S_B_it; Matrix & Piola_kirchhoff_matrix = *piola_it; SetCauchyStressMatrix< dim >(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.mul < true, false > (B, S); Bt_S_B.mul < false, false > (Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array (nb_element, bt_s_b_size * bt_s_b_size, "K_e"); model->getFEM().integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); const Array & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & strain_vect = strain(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array (nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array (nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_iterator< Matrix > shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::iterator< Matrix > shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array (nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::iterator< Matrix > Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim * nb_nodes_per_element); Array::iterator< Matrix > grad_u_it = strain_vect.begin(dim, dim); Array::iterator< Matrix > D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::iterator< Matrix > D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { Matrix & grad_u = *grad_u_it; Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; //transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.mul < true, false > (B, D); Bt_D_B.mul < false, false > (Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array (nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEM().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleResidual(GhostType ghost_type){ AKANTU_DEBUG_IN(); Array & residual = const_cast &> (model->getResidual()); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(dim, ghost_type); Mesh::type_iterator last_type = mesh.lastType(dim, ghost_type); for (; it != last_type; ++it) { const Array & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEM::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array::iterator< Matrix > shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); //Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); //Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; Array * bt_s = new Array (nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); Array::iterator< Matrix > grad_u_it = this->strain(*it, ghost_type).begin(dim, dim); Array::iterator< Matrix > grad_u_end = this->strain(*it, ghost_type).end(dim, dim); Array::iterator< Matrix > stress_it = this->piola_kirchhoff_stress(*it, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::iterator< Matrix > bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { Matrix & grad_u = *grad_u_it; Matrix & r_it = *bt_s_it; Matrix & S_it = *stress_it; SetCauchyStressArray (S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.mul < true, false > (B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array (nb_element, bt_s_size, "r_e"); model->getFEM().integrate(*bt_s, *r_e, bt_s_size, *it, ghost_type, elem_filter); delete bt_s; model->getFEM().assembleArray(*r_e, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { switch(spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(*it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & strain_vect = strain(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); Array & disp = model->getDisplacement(); model->getFEM().gradientOnQuadraturePoints(disp, strain_vect, dim, type, ghost_type, elem_filter); UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); Array * tangent_moduli_tensors = new Array(nb_element*nb_quadrature_points, tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); tangent_moduli_tensors->clear(); computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); Array * shapesd_filtered = new Array(0, dim* nb_nodes_per_element, "filtered shapesd"); FEM::filterElementalData(model->getFEM().getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array filtered_u(nb_element, nb_nodes_per_element * spatial_dimension); FEM::extractNodalToElementField(model->getFEM().getMesh(), disp, filtered_u, type, ghost_type, elem_filter); /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ Array::iterator< Matrix > shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::iterator< Matrix > D_it = tangent_moduli_tensors->begin(tangent_moduli_size, tangent_moduli_size); Array::iterator< Matrix > sigma_it = stress(type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Vector > u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); Matrix B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element); Vector Bu(tangent_moduli_size); Vector DBu(tangent_moduli_size); for (UInt e = 0; e < nb_element; ++e, ++u_it) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { Vector & u = *u_it; Matrix & sigma = *sigma_it; Matrix & D = *D_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bu.mul(B, u); DBu.mul(D, Bu); // Voigt notation to full symmetric tensor for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i); if(dim == 2) { sigma(0,1) = sigma(1,0) = DBu(2); } else if(dim == 3) { sigma(1,2) = sigma(2,1) = DBu(3); sigma(0,2) = sigma(2,0) = DBu(4); sigma(0,1) = sigma(1,0) = DBu(5); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Real * epot = potential_energy(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computePotentialEnergyOnQuad(grad_u, sigma, *epot); epot++; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { if(!potential_energy.exists(*it, _not_ghost)) { UInt nb_element = element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points){ Array::iterator< Matrix > strain_it = this->strain(type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > strain_end = this->strain(type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > stress_it = this->stress(type).begin(spatial_dimension, spatial_dimension); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type); strain_it += index*nb_quadrature_points; strain_end += (index+1)*nb_quadrature_points; stress_it += index*nb_quadrature_points; Real * epot_quad = epot_on_quad_points.storage(); for(;strain_it != strain_end; ++strain_it, ++stress_it, ++epot_quad) { Matrix & grad_u = *strain_it; Matrix & sigma = *stress_it; computePotentialEnergyOnQuad(grad_u, sigma, *epot_quad); } } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElement(); /// integrate the potential energy for each type of elements Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += model->getFEM().integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(model->getFEM().getNbQuadraturePoints(type)); computePotentialEnergyByElement(type,index,epot_on_quad_points); epot = model->getFEM().integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if(energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); Array nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += model->getDisplacement(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_tot_quad = model->getFEM().getNbQuadraturePoints(*it, ghost_type) * nb_element; Array & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); model->getFEM().interpolateOnQuadraturePoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation(const ByElementTypeReal & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); ByElementTypeArray quadrature_points_coordinates("quadrature_points_coordinates_for_interpolation", getID()); mesh.initByElementTypeArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; computeQuadraturePointsCoordinates(quadrature_points_coordinates, ghost_type); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array & interp_points_coord = interpolation_points_coordinates(type, ghost_type); UInt nb_interpolation_points_per_elem = interp_points_coord.getSize() / nb_element; AKANTU_DEBUG_ASSERT(interp_points_coord.getSize() % nb_element == 0, "Number of interpolation points is wrong"); #define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ initElementalFieldInterpolation(quadrature_points_coordinates(type, ghost_type), \ interp_points_coord, \ nb_interpolation_points_per_elem, \ ghost_type) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD); #undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const UInt nb_interpolation_points_per_elem, const GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(ghost_type); Array & elem_fil = element_filter(type, ghost_type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type, ghost_type); if(!interpolation_inverse_coordinates.exists(type, ghost_type)) interpolation_inverse_coordinates.alloc(nb_element, size_inverse_coords*size_inverse_coords, type, ghost_type); if(!interpolation_points_matrices.exists(type, ghost_type)) interpolation_points_matrices.alloc(nb_element, nb_interpolation_points_per_elem * size_inverse_coords, type, ghost_type); Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); Array & interp_points_mat = interpolation_points_matrices(type, ghost_type); Matrix quad_coord_matrix(size_inverse_coords, size_inverse_coords); Array::const_iterator< Matrix > quad_coords_it = quad_coordinates.begin_reinterpret(spatial_dimension, nb_quad_per_element, nb_element); Array::const_iterator< Matrix > points_coords_begin = interpolation_points_coordinates.begin_reinterpret(spatial_dimension, nb_interpolation_points_per_elem, interpolation_points_coordinates.getSize() / nb_interpolation_points_per_elem); Array::iterator< Matrix > inv_quad_coord_it = interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); Array::iterator< Matrix > inv_points_mat_it = interp_points_mat.begin(nb_interpolation_points_per_elem, size_inverse_coords); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++inv_quad_coord_it, ++inv_points_mat_it, ++quad_coords_it) { /// matrix containing the quadrature points coordinates const Matrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(quad_coords, quad_coord_matrix); /// invert the interpolation matrix inv_quad_coord_matrix.inverse(quad_coord_matrix); /// matrix containing the interpolation points coordinates const Matrix & points_coords = points_coords_begin[elem_fil(el)]; /// matrix to store the interpolation points coordinates /// compatible with these functions Matrix & inv_points_coord_matrix = *inv_points_mat_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(points_coords, inv_points_coord_matrix); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ByElementTypeReal & result, const GhostType ghost_type) { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; Array & res = result(type, ghost_type); #define INTERPOLATE_ELEMENTAL_FIELD(type) \ interpolateElementalField(stress(type, ghost_type), \ res, \ ghost_type) AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD); #undef INTERPOLATE_ELEMENTAL_FIELD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::interpolateElementalField(const Array & field, Array & result, const GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & elem_fil = element_filter(type, ghost_type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type, ghost_type); UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(ghost_type); Matrix coefficients(nb_quad_per_element, field.getNbComponent()); const Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); const Array & interp_points_coord = interpolation_points_matrices(type, ghost_type); UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / size_inverse_coords; Array::const_iterator< Matrix > field_it = field.begin_reinterpret(field.getNbComponent(), nb_quad_per_element, nb_element); Array::const_iterator< Matrix > interpolation_points_coordinates_it = interp_points_coord.begin(nb_interpolation_points_per_elem, size_inverse_coords); Array::iterator< Matrix > result_begin = result.begin_reinterpret(field.getNbComponent(), nb_interpolation_points_per_elem, result.getSize() / nb_interpolation_points_per_elem); Array::const_iterator< Matrix > inv_quad_coord_it = interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** * matrix containing the inversion of the quadrature points' * coordinates */ const Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// matrix containing the points' coordinates const Matrix & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the result result_begin[elem_fil(el)].mul(coefficients, coord); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const Array & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" < & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ const ByElementTypeArray & Material::getInternal(const ID & int_id) const { std::map *>::const_iterator it = internal_vectors_real.find(getID() + ":" + int_id); if(it == internal_vectors_real.end()) { AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ ByElementTypeArray & Material::getInternal(const ID & int_id) { std::map *>::iterator it = internal_vectors_real.find(getID() + ":" + int_id); if(it == internal_vectors_real.end()) { AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } +/* -------------------------------------------------------------------------- */ +void Material::addElements(const Array & elements_to_add) { + AKANTU_DEBUG_IN(); + UInt mat_id = model->getInternalIndexFromID(getID()); + Array::const_iterator el_begin = elements_to_add.begin(); + Array::const_iterator el_end = elements_to_add.end(); + for(;el_begin != el_end; ++el_begin) { + const Element & element = *el_begin; + Array & element_index_material = model->getElementIndexByMaterial(element.type, element.ghost_type); + UInt index = this->addElement(element.type, element.element, element.ghost_type); + element_index_material(element.element, 0) = mat_id; + element_index_material(element.element, 1) = index; + } + + this->resizeInternals(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Material::removeElements(const Array & elements_to_remove) { + AKANTU_DEBUG_IN(); + + Array::const_iterator el_begin = elements_to_remove.begin(); + Array::const_iterator el_end = elements_to_remove.end(); + + ByElementTypeUInt material_local_new_numbering("remove mat filter elem", getID()); + + Element element; + for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { + GhostType ghost_type = *gt; + element.ghost_type = ghost_type; + + ByElementTypeArray::type_iterator it = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); + ByElementTypeArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); + + for(; it != end; ++it) { + ElementType type = *it; + element.type = type; + + Array & elem_filter = this->element_filter(type, ghost_type); + Array & element_index_material = model->getElementIndexByMaterial(type, ghost_type); + + if(!material_local_new_numbering.exists(type, ghost_type)) + material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, ghost_type); + Array & mat_renumbering = material_local_new_numbering(type, ghost_type); + + UInt nb_element = elem_filter.getSize(); + Array elem_filter_tmp; + + UInt new_id = 0; + for (UInt el = 0; el < nb_element; ++el) { + element.element = elem_filter(el); + + if(std::find(el_begin, el_end, element) == el_end) { + elem_filter_tmp.push_back(element.element); + + mat_renumbering(el) = new_id; + element_index_material(element.element, 1) = new_id; + ++new_id; + } else { + mat_renumbering(el) = UInt(-1); + } + } + + elem_filter.resize(elem_filter_tmp.getSize()); + elem_filter.copy(elem_filter_tmp); + } + } + + for (std::map *>::iterator it = internal_vectors_real.begin(); + it != internal_vectors_real.end(); + ++it) it->second->removeQuadraturePoints(material_local_new_numbering); + + for (std::map *>::iterator it = internal_vectors_uint.begin(); + it != internal_vectors_uint.end(); + ++it) it->second->removeQuadraturePoints(material_local_new_numbering); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Material::resizeInternals() { + AKANTU_DEBUG_IN(); + for (std::map *>::iterator it = internal_vectors_real.begin(); + it != internal_vectors_real.end(); + ++it) it->second->resize(); + + for (std::map *>::iterator it = internal_vectors_uint.begin(); + it != internal_vectors_uint.end(); + ++it) it->second->resize(); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void Material::onElementsAdded(__attribute__((unused)) const Array & element_list, + __attribute__((unused)) const NewElementsEvent & event) { + this->resizeInternals(); +} + +/* -------------------------------------------------------------------------- */ +void Material::onElementsRemoved(const Array & element_list, + const ByElementTypeUInt & new_numbering, + __attribute__((unused)) const RemovedElementsEvent & event) { + UInt my_num = model->getInternalIndexFromID(getID()); + + ByElementTypeUInt material_local_new_numbering("remove mat filter elem", getID()); + + Array::const_iterator el_begin = element_list.begin(); + Array::const_iterator el_end = element_list.end(); + + for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { + GhostType gt = *g; + + ByElementTypeArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); + ByElementTypeArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, _ek_not_defined); + for (; it != end; ++it) { + ElementType type = *it; + if(element_filter.exists(type, gt)){ + Array & elem_filter = element_filter(type, gt); + + Array & element_index_material = model->getElementIndexByMaterial(type, gt); + element_index_material.resize(model->getFEM().getMesh().getNbElement(type, gt)); // all materials will resize of the same size... + + if(!material_local_new_numbering.exists(type, gt)) + material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); + + Array & mat_renumbering = material_local_new_numbering(type, gt); + const Array & renumbering = new_numbering(type, gt); + Array elem_filter_tmp; + UInt ni = 0; + Element el; + el.type = type; + el.ghost_type = gt; + for (UInt i = 0; i < elem_filter.getSize(); ++i) { + el.element = elem_filter(i); + if(std::find(el_begin, el_end, el) == el_end) { + UInt new_el = renumbering(el.element); + AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); + elem_filter_tmp.push_back(new_el); + mat_renumbering(i) = ni; + element_index_material(new_el, 0) = my_num; + element_index_material(new_el, 1) = ni; + ++ni; + } else { + mat_renumbering(i) = UInt(-1); + } + } + + elem_filter.resize(elem_filter_tmp.getSize()); + elem_filter.copy(elem_filter); + } + } + } + + for (std::map *>::iterator it = internal_vectors_real.begin(); + it != internal_vectors_real.end(); + ++it) it->second->removeQuadraturePoints(material_local_new_numbering); + + for (std::map *>::iterator it = internal_vectors_uint.begin(); + it != internal_vectors_uint.end(); + ++it) it->second->removeQuadraturePoints(material_local_new_numbering); +} + /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(":") + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 3b932b36f..7454c3c80 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,570 +1,582 @@ /** * @file material.hh * * @author Marco Vocialta * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_voigthelper.hh" #include "parser.hh" #include "parsable.hh" #include "data_accessor.hh" #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleResidual(GhostType ghost_type); /// save the stress in the precious_stress if needed virtual void savePreviousState(GhostType ghost_type); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(__attribute__((unused)) Real h, __attribute__((unused)) const Element & element = ElementNull) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); + /// add many elements at once + void addElements(const Array & elements_to_add); + + /// remove many element at once + void removeElements(const Array & elements_to_remove); + /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ void interpolateStress(ByElementTypeReal & result, const GhostType ghost_type = _not_ghost); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ void initElementalFieldInterpolation(const ByElementTypeReal & interpolation_points_coordinates); protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// assemble the residual template void assembleResidual(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation template void computeCauchyStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost); template inline void computeCauchyStressOnQuad(const Matrix & F, const Matrix & S, Matrix & cauchy); template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type); /// write the stress tensor in the Voigt notation. template inline void SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect); inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// Size of the Stress matrix for the case of finite deformation see: Bathe et al, IJNME, Vol 9, 353-386, 1975 inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const; /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, 1975 template inline void SetCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix); /// compute the potential energy by element void computePotentialEnergyByElement(); + /// resize the intenals arrays + void resizeInternals(); + public: /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const; protected: /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Array & field, Array & result, const GhostType ghost_type); /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const UInt nb_interpolation_points_per_elem, const GhostType ghost_type); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic linear elements inline void buildElementalFieldInterpolationCoodinatesLinear(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic quadratic elements inline void buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix & coordinates, Matrix & coordMatrix); /// get the size of the coordiante matrix used in the interpolation template inline UInt getSizeElementalFieldInterpolationCoodinates(GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for a quadrature point virtual inline void computePotentialEnergyOnQuad(Matrix & grad_u, Matrix & sigma, Real & epot); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points); public: /* ------------------------------------------------------------------------ */ template inline void gradUToF(const Matrix & grad_u, Matrix & F); inline void rightCauchy(const Matrix & F, Matrix & C); inline void leftCauchy (const Matrix & F, Matrix & B); template inline void deformationJacobian(const Matrix & F, Real & J); inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ByElementTypeArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ByElementTypeArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); public: /* ------------------------------------------------------------------------ */ - virtual inline void onElementsAdded(const Array & element_list, - const NewElementsEvent & event); + virtual void onElementsAdded(const Array & element_list, + const NewElementsEvent & event); - virtual inline void onElementsRemoved(const Array & element_list, - const ByElementTypeUInt & new_numbering, - const RemovedElementsEvent & event); + virtual void onElementsRemoved(const Array & element_list, + const ByElementTypeUInt & new_numbering, + const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, Memory::getID(), const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); /// return the potential energy for the subset of elements contained by the material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained by the material virtual Real getEnergy(std::string energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(Strain, strain, const ByElementTypeReal &); AKANTU_GET_MACRO(Stress, stress, const ByElementTypeReal &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ByElementTypeUInt &); bool isNonLocal() const { return is_non_local; } const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); const ByElementTypeArray & getInternal(const ID & id) const; ByElementTypeArray & getInternal(const ID & id); inline bool isInternal(const ID & id) const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); template inline const T & getParam(const ID & param) const; protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; protected: /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ByElementTypeArray element_filter; /// stresses arrays ordered by element types InternalField stress; + /// strains arrays ordered by element types + InternalField pre_strain; + /// strains arrays ordered by element types InternalField strain; /// stress increment arrays ordered by element types (Finite deformation) InternalField delta_stress; /// strain increment arrays ordered by element types (Finite deformation) InternalField delta_strain; /// inelastic strain arrays ordered by element types (inelastic deformation) InternalField inelas_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types (Finite deformation) InternalField piola_kirchhoff_stress; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// previous stresses InternalField previous_stress; /// tell if the material need the previous strain state bool use_previous_strain; /// previous strain InternalField previous_strain; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ this->stress(el_type, \ ghost_type).resize(this->strain(el_type, \ ghost_type).getSize()); \ \ Array::iterator< Matrix > stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ if(this->isFiniteDeformation()) \ Array::iterator< Matrix > stress_it = \ this->piola_kirchhoff_stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > sigma_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ tangent_mat.resize(this->strain(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Array::iterator< Matrix > tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_FINITE_DEFORMATION_BEGIN(tangent_mat) \ Array::iterator< Matrix > strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > green_it = \ this->delta_strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Array::iterator< Matrix > sigma_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ tangent_mat.resize(this->strain(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Array::iterator< Matrix > tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++green_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *strain_it; \ Matrix & __attribute__((unused)) grad_delta_u = *green_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_FINITE_DEFORMATION_END \ } \ /* -------------------------------------------------------------------------- */ #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index 48052816d..8d2f4377f 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,496 +1,419 @@ /** * @file material_inline_impl.cc * * @author Marco Vocialta * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the inline functions of the class material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ __END_AKANTU__ #include "solid_mechanics_model.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { Array & el_filter = element_filter(type, ghost_type); el_filter.push_back(element); return el_filter.getSize()-1; } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ inline UInt Material::getCauchyStressMatrixSize(UInt dim) const { return (dim * dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const Matrix & grad_u, Matrix & F) { UInt size_F = F.rows(); AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim*dim, "The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u."); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) = grad_u(i, j); for (UInt i = 0; i < size_F; ++i) F(i, i) += 1; } /* -------------------------------------------------------------------------- */ template<> inline void Material::deformationJacobian<3>(const Matrix & F, Real & J) { J = Math::det3(F.storage()); } /* -------------------------------------------------------------------------- */ template<> inline void Material::deformationJacobian<2>(const Matrix & F, Real & J) { J = Math::det2(F.storage()); } /* -------------------------------------------------------------------------- */ template<> inline void Material::deformationJacobian<1>(const Matrix & F, Real & J) { J = *(F.storage()); } /* -------------------------------------------------------------------------- */ template inline void Material::computeCauchyStressOnQuad(const Matrix & F, const Matrix & piola, Matrix & sigma) { Real J = 1.0; deformationJacobian(F, J); Matrix F_S(dim, dim); F_S.mul(F, piola); Real constant = J ? 1. /J : 0; sigma.mul(F_S, F, constant); } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const Matrix & F, Matrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const Matrix & F, Matrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) epsilon(i, j) = 0.5*(grad_u(i, j) + grad_u(j, i)); } /* -------------------------------------------------------------------------- */ inline void Material::computePotentialEnergyOnQuad(Matrix & grad_u, Matrix & sigma, Real & epot) { epot = 0.; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) epot += sigma(i, j) * grad_u(i, j); epot *= .5; } /* ---------------------------------------------------------------------------*/ template inline void Material::SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect) { AKANTU_DEBUG_IN(); Stress_vect.clear(); //UInt cauchy_matrix_size = getCauchyStressArraySize(dim); //see Finite ekement formulations for large deformation dynamic analysis, Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau /* * 1d: [ s11 ]' * 2d: [ s11 s22 s12 ]' * 3d: [ s11 s22 s33 s23 s13 s12 ] */ for (UInt i = 0; i < dim; ++i)//diagonal terms Stress_vect(i, 0) = S_t(i, i); for (UInt i = 1; i < dim; ++i)// term s12 in 2D and terms s23 s13 in 3D Stress_vect(dim+i-1, 0) = S_t(dim-i-1, dim-1); for (UInt i = 2; i < dim; ++i)//term s13 in 3D Stress_vect(dim+i, 0) = S_t(0, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Material::SetCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix) { AKANTU_DEBUG_IN(); Stress_matrix.clear(); /// see Finite ekement formulations for large deformation dynamic analysis, Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau for (UInt i = 0; i < dim; ++i) { for (UInt m = 0; m < dim; ++m) { for (UInt n = 0; n < dim; ++n) { Stress_matrix(i * dim + m, i * dim + n) = S_t(m, n); } } } //other terms from the diagonal /*for (UInt i = 0; i < 3 - dim; ++i) { Stress_matrix(dim * dim + i, dim * dim + i) = S_t(dim + i, dim + i); }*/ AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const Matrix & coordinates, __attribute__((unused)) Matrix & coordMatrix) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ inline void Material::buildElementalFieldInterpolationCoodinatesLinear(const Matrix & coordinates, Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.cols(); ++i) coordMatrix(i, 0) = 1; } /* -------------------------------------------------------------------------- */ inline void Material::buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix & coordinates, Matrix & coordMatrix) { UInt nb_quadrature_points = coordMatrix.cols(); for (UInt i = 0; i < coordinates.cols(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(j-1, i); } } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_segment_2>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_segment_3>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_tetrahedron_4>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_tetrahedron_10>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix); } /** * @todo Write a more efficient interpolation for quadrangles by * dropping unnecessary quadrature points * */ /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_4>(const Matrix & coordinates, Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.cols(); ++i) { Real x = coordinates(0, i); Real y = coordinates(1, i); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; } } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_8>(const Matrix & coordinates, Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.cols(); ++i) { UInt j = 0; Real x = coordinates(0, i); Real y = coordinates(1, i); for (UInt e = 0; e <= 2; ++e) { for (UInt n = 0; n <= 2; ++n) { coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n); ++j; } } } } /* -------------------------------------------------------------------------- */ template inline UInt Material::getSizeElementalFieldInterpolationCoodinates(GhostType ghost_type) { return model->getFEM().getNbQuadraturePoints(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataForElements(const Array & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension * spatial_dimension * sizeof(Real) * this->getModel().getNbQuadraturePoints(elements); } return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { if(this->isFiniteDeformation()) { packElementDataHelper(piola_kirchhoff_stress, buffer, elements); packElementDataHelper(strain, buffer, elements); } packElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline void Material::unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if(tag == _gst_smm_stress) { if(this->isFiniteDeformation()) { unpackElementDataHelper(piola_kirchhoff_stress, buffer, elements); unpackElementDataHelper(strain, buffer, elements); } unpackElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ template inline const T & Material::getParam(const ID & param) const { try { return get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { set(param, value); } catch(...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template inline void Material::packElementDataHelper(const ByElementTypeArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) const { DataAccessor::packElementalDataHelper(data_to_pack, buffer, elements, true, model->getFEM(fem_id)); } /* -------------------------------------------------------------------------- */ template inline void Material::unpackElementDataHelper(ByElementTypeArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) { DataAccessor::unpackElementalDataHelper(data_to_unpack, buffer, elements, true, model->getFEM(fem_id)); } - -/* -------------------------------------------------------------------------- */ -inline void Material::onElementsAdded(__attribute__((unused)) const Array & element_list, - __attribute__((unused)) const NewElementsEvent & event) { - for (std::map *>::iterator it = internal_vectors_real.begin(); - it != internal_vectors_real.end(); - ++it) it->second->resize(); - - for (std::map *>::iterator it = internal_vectors_uint.begin(); - it != internal_vectors_uint.end(); - ++it) it->second->resize(); -} - -/* -------------------------------------------------------------------------- */ -inline void Material::onElementsRemoved(const Array & element_list, - const ByElementTypeUInt & new_numbering, - __attribute__((unused)) const RemovedElementsEvent & event) { - UInt my_num = model->getInternalIndexFromID(getID()); - - ByElementTypeUInt material_local_new_numbering("remove mat filter elem", getID()); - - Array::const_iterator el_begin = element_list.begin(); - Array::const_iterator el_end = element_list.end(); - - for(UInt g = _not_ghost; g <= _ghost; ++g) { - GhostType gt = (GhostType) g; - ByElementTypeArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); - ByElementTypeArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, _ek_not_defined); - for (; it != end; ++it) { - ElementType type = *it; - if(element_filter.exists(type, gt)){ - Array & elem_filter = element_filter(type, gt); - - Array & element_index_material = model->getElementIndexByMaterial(type, gt); - element_index_material.resize(model->getFEM().getMesh().getNbElement(type, gt)); // all materials will resize of the same size... - - if(!material_local_new_numbering.exists(type, gt)) - material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); - - Array & mat_renumbering = material_local_new_numbering(type, gt); - const Array & renumbering = new_numbering(type, gt); - Array elem_filter_tmp; - UInt ni = 0; - Element el; - el.type = type; - el.ghost_type = gt; - for (UInt i = 0; i < elem_filter.getSize(); ++i) { - el.element = elem_filter(i); - if(std::find(el_begin, el_end, el) == el_end) { - UInt new_el = renumbering(el.element); - AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); - elem_filter_tmp.push_back(new_el); - mat_renumbering(i) = ni; - element_index_material(new_el, 0) = my_num; - element_index_material(new_el, 1) = ni; - ++ni; - } else { - mat_renumbering(i) = UInt(-1); - } - } - - elem_filter.resize(elem_filter_tmp.getSize()); - elem_filter.copy(elem_filter); - } - } - } - - for (std::map *>::iterator it = internal_vectors_real.begin(); - it != internal_vectors_real.end(); - ++it) it->second->removeQuadraturePoints(material_local_new_numbering); - - for (std::map *>::iterator it = internal_vectors_uint.begin(); - it != internal_vectors_uint.end(); - ++it) it->second->removeQuadraturePoints(material_local_new_numbering); -} - -/* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template<> inline void Material::registerInternal(InternalField & vect) { internal_vectors_real[vect.getID()] = &vect; } template<> inline void Material::registerInternal(InternalField & vect) { internal_vectors_uint[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template<> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_real.erase(vect.getID()); } template<> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_uint.erase(vect.getID()); } /* -------------------------------------------------------------------------- */ inline bool Material::isInternal(const ID & id) const { return internal_vectors_real.find(this->getID()+":"+id) != internal_vectors_real.end(); } diff --git a/src/model/solid_mechanics/material_selector_tmpl.hh b/src/model/solid_mechanics/material_selector_tmpl.hh index 1e3047315..260342f3f 100644 --- a/src/model/solid_mechanics/material_selector_tmpl.hh +++ b/src/model/solid_mechanics/material_selector_tmpl.hh @@ -1,87 +1,87 @@ /** * @file material_selector_tmpl.hh * * @author Nicolas Richart * * @date Tue Aug 13 11:06:50 2013 * * @brief Implementation of the template MaterialSelector * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_SELECTOR_TMPL_HH__ #define __AKANTU_MATERIAL_SELECTOR_TMPL_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template<> class MeshDataMaterialSelector : public MaterialSelector { public: MeshDataMaterialSelector(const std::string & name, const SolidMechanicsModel & model) : names(model.getMesh().getData(name)), model(model) { } UInt operator() (const Element & element) { try { DebugLevel dbl = debug::getDebugLevel(); debug::setDebugLevel(dblError); std::string material_name = names(element.type, element.ghost_type)(element.element); debug::setDebugLevel(dbl); return model.getMaterialIndex(material_name); } catch (...) { return MaterialSelector::operator()(element); } } private: const ByElementTypeArray & names; const SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ template<> class MeshDataMaterialSelector : public MaterialSelector { public: MeshDataMaterialSelector(const std::string & name, const SolidMechanicsModel & model, UInt first_index = 1) : indexes(model.getMesh().getData(name)), model(model), first_index(first_index) { } UInt operator() (const Element & element) { try { DebugLevel dbl = debug::getDebugLevel(); debug::setDebugLevel(dblError); UInt mat = indexes(element.type, element.ghost_type)(element.element) - first_index; debug::setDebugLevel(dbl); return mat; } catch (...) { return MaterialSelector::operator()(element); } } -private: +protected: const ByElementTypeArray indexes; const SolidMechanicsModel & model; UInt first_index; }; __END_AKANTU__ #endif /* __AKANTU_MATERIAL_SELECTOR_TMPL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage.hh b/src/model/solid_mechanics/materials/material_damage/material_damage.hh index 14f5801bf..03032d8bd 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_damage.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_damage.hh @@ -1,103 +1,119 @@ /** * @file material_damage.hh * * @author Marion Estelle Chambart * @author Nicolas Richart * * @date Tue Mar 15 16:06:20 2011 * * @brief Material isotropic elastic * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_HH__ #define __AKANTU_MATERIAL_DAMAGE_HH__ __BEGIN_AKANTU__ template class Parent = MaterialElastic> class MaterialDamage : public Parent { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDamage(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialDamage() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual void computeAllStresses(GhostType ghost_type); + /// constitutive law for all element of a type + virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); + + /// compute the tangent stiffness matrix for an element type + virtual void computeTangentModuli(const ElementType & el_type, + Array & tangent_matrix, + GhostType ghost_type = _not_ghost); + protected: /// update the dissipated energy, must be called after the stress have been computed void updateDissipatedEnergy(GhostType ghost_type); + /// constitutive law for a given quadrature point + inline void computeStressOnQuad(const Matrix & grad_u, + Matrix & sigma, + Real & dam); + + /// compute the tangent stiffness matrix for a given quadrature point + inline void computeTangentModuliOnQuad(Matrix & tangent, Real & dam); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// give the dissipated energy for the time step Real getDissipatedEnergy() const; virtual Real getEnergy(std::string type); virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) { return Parent::getEnergy(energy_id, type, index); }; AKANTU_GET_MACRO_NOT_CONST(Damage, damage, ByElementTypeReal &); AKANTU_GET_MACRO(Damage, damage, const ByElementTypeReal &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// damage internal variable InternalField damage; /// dissipated energy InternalField dissipated_energy; /// contain the current value of @f$ \int_0^{\epsilon}\sigma(\omega)d\omega @f$ the dissipated energy InternalField int_sigma; + }; __END_AKANTU__ #include "material_damage_tmpl.hh" #endif /* __AKANTU_MATERIAL_DAMAGE_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh b/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh index de0fbc22d..3f6f2e675 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh @@ -1,154 +1,210 @@ /** * @file material_damage.cc * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date Tue Mar 15 16:06:20 2011 * * @brief Specialization of the material class for the damage material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_damage.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class Parent> MaterialDamage::MaterialDamage(SolidMechanicsModel & model, const ID & id) : Material(model, id), Parent(model, id), damage("damage", *this), dissipated_energy("damage dissipated energy", *this), int_sigma("integral of sigma", *this) { AKANTU_DEBUG_IN(); this->is_non_local = false; this->use_previous_stress = true; this->use_previous_strain = true; this->damage .initialize(1); this->dissipated_energy.initialize(1); this->int_sigma .initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialDamage::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the dissipated energy in each element by a trapezoidal approximation * of * @f$ Ed = \int_0^{\epsilon}\sigma(\omega)d\omega - \frac{1}{2}\sigma:\epsilon@f$ */ template class Parent> void MaterialDamage::updateDissipatedEnergy(GhostType ghost_type) { // compute the dissipated energy per element const Mesh & mesh = this->model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { ElementType el_type = *it; Array::iterator< Matrix > sigma = this->stress(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > sigma_p = this->previous_stress(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > epsilon = this->strain(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > epsilon_p = this->previous_strain(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator ints = int_sigma(el_type, ghost_type).begin(); Array::iterator ed = dissipated_energy(el_type, ghost_type).begin(); Array::iterator ed_end = dissipated_energy(el_type, ghost_type).end(); for (; ed != ed_end; ++ed, ++ints, ++epsilon, ++sigma, ++epsilon_p, ++sigma_p) { Real epot = 0.; Real dint = 0.; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { epot += (*sigma)(i,j) * (*epsilon)(i,j); /// \f$ epot = .5 \sigma : \epsilon \f$ dint += .5 * ((*sigma_p)(i,j) + (*sigma)(i,j)) * ((*epsilon)(i,j) - (*epsilon_p)(i,j)); /// \f$ \frac{.5 \sigma(\epsilon(t-h)) + \sigma(\epsilon(t))}{\epsilon(t) - \epsilon(t-h)} \f$ (*epsilon_p)(i,j) = (*epsilon)(i,j); (*sigma_p)(i,j) = (*sigma)(i,j); } } epot *= .5; *ints += dint; *ed = *ints - epot; } } } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialDamage::computeAllStresses(GhostType ghost_type) { Material::computeAllStresses(ghost_type); if(!this->is_non_local) this->updateDissipatedEnergy(ghost_type); } +/* -------------------------------------------------------------------------- */ +template class Parent> +void MaterialDamage::computeStress(ElementType el_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + + Real * dam = this->damage(el_type, ghost_type).storage(); + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + computeStressOnQuad(grad_u, sigma, *dam); + + ++dam; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template class Parent> +void MaterialDamage::computeTangentModuli(const ElementType & el_type, + Array & tangent_matrix, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + Real * dam = this->damage(el_type, ghost_type).storage(); + + MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); + computeTangentModuliOnQuad(tangent, *dam); + + ++dam; + + MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template class Parent> +void MaterialDamage::computeStressOnQuad(const Matrix & grad_u, + Matrix & sigma, + + Real & dam) { + Parent::computeStressOnQuad(grad_u, sigma); + sigma *= (1-dam); + +} + +/* -------------------------------------------------------------------------- */ +template class Parent> +void MaterialDamage::computeTangentModuliOnQuad(Matrix & tangent, + Real & dam) { + Parent::computeTangentModuliOnQuad(tangent); + tangent *= (1-dam); +} + /* -------------------------------------------------------------------------- */ template class Parent> Real MaterialDamage::getDissipatedEnergy() const { AKANTU_DEBUG_IN(); Real de = 0.; const Mesh & mesh = this->model->getFEM().getMesh(); /// integrate the dissipated energy for each type of elements Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost); for(; it != end; ++it) { de += this->model->getFEM().integrate(dissipated_energy(*it, _not_ghost), *it, _not_ghost, this->element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return de; } /* -------------------------------------------------------------------------- */ template class Parent> Real MaterialDamage::getEnergy(std::string type) { if(type == "dissipated") return getDissipatedEnergy(); else return Parent::getEnergy(type); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_damage/material_mazars.cc b/src/model/solid_mechanics/materials/material_damage/material_mazars.cc index df405cbd5..b651d9e91 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_mazars.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_mazars.cc @@ -1,82 +1,82 @@ /** * @file material_mazars.cc * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date Wed Apr 06 10:09:38 2011 * * @brief Specialization of the material class for the damage material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_mazars.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialMazars::MaterialMazars(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamage(model, id), K0("K0", *this), damage_in_compute_stress(true) { AKANTU_DEBUG_IN(); this->registerParam("K0" , K0 , _pat_parsable, "K0"); this->registerParam("At" , At , 0.8 , _pat_parsable, "At"); this->registerParam("Ac" , Ac , 1.4 , _pat_parsable, "Ac"); this->registerParam("Bc" , Bc , 1900. , _pat_parsable, "Bc"); this->registerParam("Bt" , Bt , 12000., _pat_parsable, "Bt"); this->registerParam("beta", beta, 1.06 , _pat_parsable, "beta"); this->K0.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazars::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Real Ehat = 0; computeStressOnQuad(grad_u, sigma, *dam, Ehat); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialMazars); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_elastic.hh b/src/model/solid_mechanics/materials/material_elastic.hh index a88779574..f89a6a219 100644 --- a/src/model/solid_mechanics/materials/material_elastic.hh +++ b/src/model/solid_mechanics/materials/material_elastic.hh @@ -1,127 +1,127 @@ /** * @file material_elastic.hh * * @author Nicolas Richart * @author Lucas Frerot * * @date Wed Aug 04 10:58:42 2010 * * @brief Material isotropic elastic * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material_thermal.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ELASTIC_HH__ #define __AKANTU_MATERIAL_ELASTIC_HH__ __BEGIN_AKANTU__ /** * Material elastic isotropic * * parameters in the material files : * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - Plane_Stress : if 0: plane strain, else: plane stress (default: 0) */ template class MaterialElastic : public MaterialThermal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialElastic(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialElastic() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type - void computeTangentModuli(const ElementType & el_type, + virtual void computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() const; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() const; protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(const Matrix & grad_u, Matrix & sigma, const Real sigma_th = 0); /// compute the tangent stiffness matrix for an element - void computeTangentModuliOnQuad(Matrix & tangent); + inline void computeTangentModuliOnQuad(Matrix & tangent); /// recompute the lame coefficient if E or nu changes virtual void updateInternalParameters(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the stable time step inline Real getStableTimeStep(Real h, const Element & element); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// First Lamé coefficient Real lambda; /// Second Lamé coefficient (shear modulus) Real mu; /// Bulk modulus Real kpa; /// Plane stress or plane strain bool plane_stress; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_elastic_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_ELASTIC_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_extra_includes.hh b/src/model/solid_mechanics/materials/material_extra_includes.hh index 5227b4e4a..059e6981e 100644 --- a/src/model/solid_mechanics/materials/material_extra_includes.hh +++ b/src/model/solid_mechanics/materials/material_extra_includes.hh @@ -1,68 +1,70 @@ /** * @file material_extra_includes.hh * * @author Nicolas Richart * * @date Wed Oct 31 16:24:42 2012 * * @brief Extra list of materials * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ // visco-elastic materials #include "material_stiffness_proportional.hh" #include "material_standard_linear_solid_deviatoric.hh" // elastic materials #include "material_neohookean.hh" #include "material_elastic_orthotropic.hh" #include "material_elastic_linear_anisotropic.hh" // damage materials +#include "material_damage.hh" #include "material_marigo.hh" #include "material_mazars.hh" #include "material_damage_linear.hh" #include "material_vreepeerlings.hh" // plasticity #include "material_plastic.hh" #include "material_plasticityinc.hh" #include "material_viscoplasticity.hh" #define AKANTU_EXTRA_MATERIAL_LIST \ ((2, (damage_linear , MaterialDamageLinear ))) \ + ((2, (damage , MaterialDamage ))) \ ((2, (marigo , MaterialMarigo ))) \ ((2, (mazars , MaterialMazars ))) \ ((2, (vreepeerlings , MaterialVreePeerlings ))) \ ((2, (ve_stiffness_prop , MaterialStiffnessProportional ))) \ ((2, (sls_deviatoric , MaterialStandardLinearSolidDeviatoric))) \ ((2, (neohookean , MaterialNeohookean ))) \ ((2, (elastic_orthotropic, MaterialElasticOrthotropic ))) \ ((2, (anisotropic , MaterialElasticLinearAnisotropic ))) \ ((2, (plastic , MaterialPlastic ))) \ ((2, (plasticity , MaterialPlasticityinc ))) \ ((2, (viscoplasticity , MaterialViscoPlasticity ))) diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 41bc1aee7..8b6736853 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1891 +1,1935 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" # include "dumper_iohelper_tmpl.hh" # include "dumper_iohelper_tmpl_homogenizing_field.hh" # include "dumper_iohelper_tmpl_material_internal_field.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), Dumpable(), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), element_index_by_material("element index by material", id), material_selector(new DefaultMaterialSelector(element_index_by_material)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->registerDumper("paraview_all", id, true); this->addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; if(solver) delete solver; if(mass_matrix) delete mass_matrix; if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; if(synch_parallel) delete synch_parallel; if(is_default_material_selector) delete material_selector; AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(std::string input_file, const ModelOptions & options) { Model::initFull(input_file, options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } // initialize the materials if(input_file != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEMBoundary() { FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnControlPoints(_not_ghost); fem_boundary.computeNormalsOnControlPoints(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { AKANTU_DEBUG_IN(); //in case of switch from implicit to explicit if(!this->isExplicit()) method = analysis_method; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; // std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); element_index_by_material.alloc(nb_element, 2, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->values; Real * position_val = mesh.getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext if(need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.savePreviousState(_not_ghost); mat.savePreviousState(_ghost); mat.computeAllStresses(_not_ghost); if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant non local stresses"); synch_registry->waitEndSynchronize(_gst_mnl_for_average); AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.savePreviousState(_not_ghost); mat.savePreviousState(_ghost); mat.computeAllStresses(_not_ghost); if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } boundary_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *boundary, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & boundary, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * boundary_val = boundary.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*boundary_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; boundary_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { if(previous_displacement) increment->copy(*previous_displacement); else increment->copy(*displacement); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); } #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS if(solver) solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; if (!increment) setIncrementFlagOn(); initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*boundary); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); // updateResidualInternal(); solve(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); increment->clear(); solver->setRHS(*residual); solver->solve(*increment); Real * increment_val = increment->storage(); Real * displacement_val = displacement->storage(); bool * boundary_val = boundary->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } else { *increment_val = 0.0; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic(Array & boundary_normal, Array & EulerAngles) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); // if(method != _static) jacobian_matrix->copyContent(*stiffness_matrix); Array * residual_rotated = new Array (nb_nodes, nb_degree_of_freedom, "residual_rotated"); //stiffness_matrix->saveMatrix("stiffness_original.out"); jacobian_matrix->applyBoundaryNormal(boundary_normal, EulerAngles, *residual, (*stiffness_matrix).getA(), *residual_rotated); //jacobian_matrix->saveMatrix("stiffness_rotated_dir.out"); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual_rotated); delete residual_rotated; if (!increment) setIncrementFlagOn(); solver->solve(*increment); Matrix T(nb_degree_of_freedom, nb_degree_of_freedom); Matrix small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Matrix T_small_rhs(nb_degree_of_freedom, nb_degree_of_freedom); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool constrain_ij = false; for (UInt j = 0; j < nb_degree_of_freedom; j++) { if (boundary_normal(n, j)) { constrain_ij = true; break; } } if (constrain_ij) { if (nb_degree_of_freedom == 2) { Real Theta = EulerAngles(n, 0); T(0, 0) = cos(Theta); T(0, 1) = -sin(Theta); T(1, 1) = cos(Theta); T(1, 0) = sin(Theta); } else if (nb_degree_of_freedom == 3) { Real Theta_x = EulerAngles(n, 0); Real Theta_y = EulerAngles(n, 1); Real Theta_z = EulerAngles(n, 2); T(0, 0) = cos(Theta_y) * cos(Theta_z); T(0, 1) = -cos(Theta_y) * sin(Theta_z); T(0, 2) = sin(Theta_y); T(1, 0) = cos(Theta_x) * sin(Theta_z) + cos(Theta_z) * sin(Theta_x) * sin(Theta_y); T(1, 1) = cos(Theta_x) * cos(Theta_z) - sin(Theta_x) * sin(Theta_y) * sin(Theta_z); T(1, 2) = -cos(Theta_y) * sin(Theta_x); T(2, 0) = sin(Theta_x) * sin(Theta_z) - cos(Theta_x) * cos(Theta_z) * sin(Theta_y); T(2, 1) = cos(Theta_z) * sin(Theta_x) + cos(Theta_x) * sin(Theta_y) * sin(Theta_z); T(2, 2) = cos(Theta_x) * cos(Theta_y); } small_rhs.clear(); T_small_rhs.clear(); for (UInt j = 0; j < nb_degree_of_freedom; j++) if(!(boundary_normal(n,j)) ) small_rhs(j,j)=increment_val[j]; T_small_rhs.mul(T,small_rhs); for (UInt j = 0; j < nb_degree_of_freedom; j++){ if(!(boundary_normal(n,j))){ for (UInt k = 0; k < nb_degree_of_freedom; k++) displacement_val[k]+=T_small_rhs(k,j); } } displacement_val += nb_degree_of_freedom; boundary_val += nb_degree_of_freedom; increment_val += nb_degree_of_freedom; } else { for (UInt j = 0; j < nb_degree_of_freedom; j++, ++displacement_val, ++increment_val, ++boundary_val) { if (!(*boundary_val)) { *displacement_val += *increment_val; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * boundary_val = boundary->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; return (error < tolerance); } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); if(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->values; Real * disp_val = displacement->values; bool * boun_val = boundary->values; for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ *incr_val *= (1. - *boun_val); *disp_val += *incr_val; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateIncrement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->values; Real * disp_val = displacement->values; Real * prev_disp_val = previous_displacement->values; for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val) *incr_val = (*disp_val - *prev_disp_val); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updatePreviousDisplacement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::iterator< Vector > eibm = element_index_by_material(*it, ghost_type).begin(2); Array X(0, nb_nodes_per_element*spatial_dimension); FEM::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::iterator< Matrix > X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) { elem.element = el; Real el_size = getFEM().getElementInradius(*X_el, *it); Real el_dt = mat_val[(*eibm)(0)]->getStableTimeStep(el_size, elem); min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEM().interpolateOnQuadraturePoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::iterator< Vector > vit = vel_on_quad.begin(spatial_dimension); Array::iterator< Vector > vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[element_index_by_material(type)(index, 0)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEM().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = boundary->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index){ AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; Vector mat = element_index_by_material(type, _not_ghost).begin(2)[index]; Real energy = materials[mat(0)]->getEnergy(energy_id, type, mat(1)); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(boundary ) boundary ->resize(nb_nodes); if(previous_displacement) previous_displacement->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } if (method != _explicit_lumped_mass) { delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material UInt mat_id = 0; UInt * mat_id_vect = NULL; try { const NewMaterialElementsEvent & event_mat = dynamic_cast(event); mat_id_vect = event_mat.getMaterialList().storage(); } catch(...) { } for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(mat_id_vect) mat_id = mat_id_vect[el]; else mat_id = (*material_selector)(elem); Material & mat = *materials[mat_id]; UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); Vector id(2); id[0] = mat_id; id[1] = mat_index; if(!element_index_by_material.exists(elem.type, elem.ghost_type)) element_index_by_material.alloc(0, 2, elem.type, elem.ghost_type); element_index_by_material(elem.type, elem.ghost_type).push_back(id); } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } if(method != _explicit_lumped_mass) AKANTU_DEBUG_TO_IMPLEMENT(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event) { // MeshUtils::purifyMesh(mesh); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(boundary ) mesh.removeNodesFromArray(*boundary , new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); } +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::reassignMaterial() { + AKANTU_DEBUG_IN(); + + std::vector< Array > element_to_add (materials.size()); + std::vector< Array > element_to_remove(materials.size()); + + Element element; + for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { + GhostType ghost_type = *gt; + element.ghost_type = ghost_type; + + Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); + Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); + for(; it != end; ++it) { + ElementType type = *it; + element.type = type; + + UInt nb_element = mesh.getNbElement(type, ghost_type); + Array & el_index_by_mat = element_index_by_material(type, ghost_type); + + for (UInt el = 0; el < nb_element; ++el) { + element.element = el; + + UInt old_material = el_index_by_mat(el, 0); + UInt new_material = (*material_selector)(element); + + if(old_material != new_material) { + element_to_add [new_material].push_back(element); + element_to_remove[old_material].push_back(element); + } + } + } + } + + std::vector::iterator mat_it; + UInt mat_index = 0; + for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it, ++mat_index) { + (*mat_it)->removeElements(element_to_remove[mat_index]); + (*mat_it)->addElements (element_to_add[mat_index]); + } + + AKANTU_DEBUG_OUT(); +} /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ internalAddDumpFieldToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField<>(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else if(field_id == "element_index_by_material") { internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular)); } else { bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_id); } if (is_internal) { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::HomogenizedField (*this, field_id, spatial_dimension, _not_ghost, _ek_regular)); } else { try { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::ElementalField(mesh.getData(field_id), spatial_dimension, _not_ghost, _ek_regular)); } catch (...) {} try { internalAddDumpFieldToDumper (dumper_name, field_id, new DumperIOHelper::ElementalField(mesh.getData(field_id), spatial_dimension, _not_ghost, _ek_regular)); } catch (...) {} } } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ByElementTypeArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ group.addDumpFieldExternalToDumper(dumper_name, \ BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "increment" ) { ADD_FIELD(increment , Real); } else if(field_id == "partitions" ) { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementPartitionField(mesh, dim, _not_ghost, _ek_regular, elemental_filter)); } else if(field_id == "element_index_by_material") { group.addDumpFieldExternalToDumper(dumper_name, field_id, new DumperIOHelper::ElementalField(element_index_by_material, dim, _not_ghost, _ek_regular, elemental_filter)); } else { typedef DumperIOHelper::HomogenizedField Field; group.addDumpFieldExternalToDumper(dumper_name, field_id, new Field(*this, field_id, dim, _not_ghost, _ek_regular, elemental_filter)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupField(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); group.removeDumpFieldFromDumper(dumper_name, field_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ f->setPadding(3); \ internalAddDumpFieldToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name) { ElementGroup & group = mesh.getElementGroup(group_name); this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(), field_id, group_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name) { #ifdef AKANTU_USE_IOHELPER ElementGroup & group = mesh.getElementGroup(group_name); UInt dim = group.getDimension(); const Array * nodal_filter = &(group.getNodes()); const ByElementTypeArray * elemental_filter = &(group.getElements()); #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field, 0, 0, nodal_filter); \ f->setPadding(3); \ group.addDumpFieldExternalToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, dim, _not_ghost, _ek_regular, elemental_filter); field->setPadding(3); group.addDumpFieldExternalToDumper(dumper_name, field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER if(field_id == "stress") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else if(field_id == "strain") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); internalAddDumpFieldToDumper(dumper_name, field_id, field); } #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_index_by_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index daa03b1e9..5c7f929ae 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,671 +1,675 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "model.hh" #include "data_accessor.hh" #include "mesh.hh" #include "dumpable.hh" #include "boundary_condition.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "integration_scheme_2nd_order.hh" #include "solver.hh" #include "material_selector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class IntegrationScheme2ndOrder; class Contact; class SparseMatrix; } __BEGIN_AKANTU__ struct SolidMechanicsModelOptions : public ModelOptions { SolidMechanicsModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, bool no_init_materials = false) : - analysis_method(analysis_method), no_init_materials(no_init_materials) {} + analysis_method(analysis_method), + no_init_materials(no_init_materials) {} AnalysisMethod analysis_method; bool no_init_materials; }; extern const SolidMechanicsModelOptions default_solid_mechanics_model_options; class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler, public Dumpable, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; typedef FEMTemplate MyFEMType; SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(std::string material_file, const ModelOptions & options = default_solid_mechanics_model_options); /// initialize the fem object needed for boundary conditions void initFEMBoundary(); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// allocate all vectors void initArrays(); /// allocate all vectors void initArraysPreviousDisplacment(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(AnalysisMethod analysis_method = _explicit_lumped_mass); bool isExplicit() { return method == _explicit_lumped_mass|| method == _explicit_consistent_mass; } /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme virtual void updateResidual(bool need_initialize = true); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); void updateIncrement(); void updatePreviousDisplacement(); /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix void solveLumped(Array & x, const Array & A, const Array & b, const Array & boundary, Real alpha); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); public: void solveStep(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false, SolverOptions & solver_options = _solver_no_options); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); public: template bool solveStep(Real tolerance, UInt max_iteration = 100); public: /// solve @f[ A\delta u = f_ext - f_int @f] in displacement void solveDynamic(); /// solve Ku = f void solveStatic(); /// solve Ku = f void solveStatic(Array & boundary_normal, Array & EulerAngles); /// test if the system is converged template bool testConvergence(Real tolerance, Real & error); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance) __attribute__((deprecated)); bool testConvergenceIncrement(Real tolerance, Real & error) __attribute__((deprecated)); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance) __attribute__((deprecated)); bool testConvergenceResidual(Real tolerance, Real & error) __attribute__((deprecated)); /// create and return the velocity damping matrix SparseMatrix & initVelocityDampingMatrix(); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); protected: /// finish the computation of residual to solve in increment void updateResidualInternal(); /// compute the support reaction and store it in force void updateSupportReaction(); public: //protected: Daniel changed it just for a test /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template void solve(Array & increment, Real block_val = 1.); /* ------------------------------------------------------------------------ */ /* Explicit/Implicit */ /* ------------------------------------------------------------------------ */ public: /// Update the stresses for the computation of the residual of the Stiffness /// matrix in the case of finite deformation void computeStresses(); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// registers all the custom materials of a given type present in the input file template void registerNewCustomMaterials(const ID & mat_type); /// register an empty material of a given type template Material & registerNewEmptyMaterial(const std::string & name); // /// Use a UIntData in the mesh to specify the material to use per element // void setMaterialIDsFromIntData(const std::string & data_name); + /// reassigns materials depending on the material selector + void reassignMaterial(); + protected: /// register a material in the dynamic database template Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline virtual void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline virtual void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); protected: inline void splitElementByMaterial(const Array & elements, Array * elements_per_mat) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded (const Array & nodes_list, const NewNodesEvent & event); virtual void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event); virtual void onElementsAdded (const Array & nodes_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupField(const std::string & field_id, const std::string & group_name); virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id); virtual void addDumpGroupFieldVector(const std::string & field_id, const std::string & group_name); virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name); virtual void addDumpFieldTensorToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step void setTimeStep(Real time_step); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition AKANTU_GET_MACRO(CurrentPosition, *current_position, const Array &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Array &); /// get the SolidMechanicsModel::residual vector, computed by /// SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, Array &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Array &); /// get the SolidMechnicsModel::incrementAcceleration vector AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Array &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); }; inline void setMaterialSelector(MaterialSelector & selector); /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be given too) Real getExternalWork(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the global jacobian matrix of the system AKANTU_GET_MACRO(GlobalJacobianMatrix, *jacobian_matrix, const SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); /// get the velocity damping matrix AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, SparseMatrix &); /// get the FEM object to integrate or interpolate on the boundary inline FEM & getFEMBoundary(const ID & name = ""); /// get integrator AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &); /// get access to the internal solver AKANTU_GET_MACRO(Solver, *solver, Solver &); /// get synchronizer AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const DistributedSynchronizer &); AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeArray &); /// vectors containing local material element index for each global element index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); protected: friend class Material; template class WeightFunction> friend class MaterialNonLocal; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement; /// lumped mass array Array * mass; /// velocities array Array * velocity; /// accelerations array Array * acceleration; /// accelerations array Array * increment_acceleration; /// forces array Array * force; /// residuals array Array * residual; /// boundaries array Array * boundary; /// array of current position used during update residual Array * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta /// t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// vectors containing local material element index for each global element index ByElementTypeUInt element_index_by_material; /// list of used materials std::vector materials; /// mapping between material name and material internal id std::map materials_names_to_id; /// class defining of to choose a material MaterialSelector * material_selector; /// define if it is the default selector or not bool is_default_material_selector; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Array * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// internal synchronizer for parallel computations DistributedSynchronizer * synch_parallel; /// tells if the material are instantiated bool are_materials_instantiated; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" #include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "solid_mechanics_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/CMakeLists.txt index d679b4959..6b319722e 100644 --- a/test/test_model/test_solid_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/CMakeLists.txt @@ -1,181 +1,196 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date Fri Sep 03 15:56:56 2010 # # @brief configuratio for SolidMechanicsModel tests # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_akantu_test(test_materials "test_materials") add_akantu_test(patch_tests "patch_tests") add_akantu_test(test_cohesive "cohesive_test" PACKAGE cohesive_element) add_akantu_test(test_contact "test_contact" PACKAGE dead_contact) #=============================================================================== add_mesh(test_solid_mechanics_model_square_mesh square.geo 2 1) add_mesh(test_solid_mechanics_model_circle_mesh1 circle.geo 2 1 OUTPUT circle1.msh) add_mesh(test_solid_mechanics_model_circle_mesh2 circle.geo 2 2 OUTPUT circle2.msh) register_test(test_solid_mechanics_model_square SOURCES test_solid_mechanics_model_square.cc DEPENDENCIES test_solid_mechanics_model_square_mesh FILES_TO_COPY material.dat test_cst_energy.pl DIRECTORIES_TO_CREATE paraview ) register_test(test_solid_mechanics_model_circle_2 SOURCES test_solid_mechanics_model_circle_2.cc DEPENDENCIES test_solid_mechanics_model_circle_mesh2 FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) #=============================================================================== add_mesh(test_bar_traction_2d_mesh1 bar.geo 2 1 OUTPUT bar1.msh) add_mesh(test_bar_traction_2d_mesh2 bar.geo 2 2 OUTPUT bar2.msh) add_mesh(test_bar_traction_2d_mesh_structured1 bar_structured.geo 2 1 OUTPUT bar_structured1.msh) register_test(test_solid_mechanics_model_bar_traction2d SOURCES test_solid_mechanics_model_bar_traction2d.cc DEPENDENCIES test_bar_traction_2d_mesh1 test_bar_traction_2d_mesh2 FILES_TO_COPY material.dat test_cst_energy.pl DIRECTORIES_TO_CREATE paraview ) register_test(test_solid_mechanics_model_bar_traction2d_structured SOURCES test_solid_mechanics_model_bar_traction2d_structured.cc DEPENDENCIES test_bar_traction_2d_mesh_structured1 FILES_TO_COPY material.dat test_cst_energy.pl DIRECTORIES_TO_CREATE paraview ) #=============================================================================== add_mesh(test_solid_mechanics_model_segment_mesh segment.geo 1 2) register_test(test_solid_mechanics_model_bar_traction2d_parallel SOURCES test_solid_mechanics_model_bar_traction2d_parallel.cc DEPENDENCIES test_bar_traction_2d_mesh2 FILES_TO_COPY material.dat test_cst_energy.pl DIRECTORIES_TO_CREATE paraview ) register_test(test_solid_mechanics_model_segment_parallel SOURCES test_solid_mechanics_model_segment_parallel.cc DEPENDENCIES test_solid_mechanics_model_segment_mesh FILES_TO_COPY material.dat test_cst_energy.pl DIRECTORIES_TO_CREATE paraview ) #=============================================================================== #register_test(test_solid_mechanics_model_bar_traction2d_mass_not_lumped # SOURCES test_solid_mechanics_model_bar_traction2d_mass_not_lumped.cc # DEPENDENCIES test_bar_traction_2d_mesh1 test_bar_traction_2d_mesh2 # FILES_TO_COPY material.dat # DIRECTORIES_TO_CREATE paraview # ) #=============================================================================== add_mesh(test_solid_mechanics_model_segment_mesh1 segment.geo 1 1 OUTPUT segment1.msh) add_mesh(test_implicit_mesh1 square_implicit.geo 2 1 OUTPUT square_implicit1.msh) add_mesh(test_implicit_mesh2 square_implicit.geo 2 2 OUTPUT square_implicit2.msh) register_test(test_solid_mechanics_model_implicit_1d SOURCES test_solid_mechanics_model_implicit_1d.cc DEPENDENCIES test_solid_mechanics_model_segment_mesh1 FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) register_test(test_solid_mechanics_model_implicit_2d SOURCES test_solid_mechanics_model_implicit_2d.cc DEPENDENCIES test_implicit_mesh1 test_implicit_mesh2 FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) #=============================================================================== add_mesh(test_implicit_beam_2d_1 beam_2d.geo 2 1 OUTPUT beam_2d_lin.msh) add_mesh(test_implicit_beam_2d_2 beam_2d.geo 2 2 OUTPUT beam_2d_quad.msh) add_mesh(test_implicit_beam_3d_2 beam_3d.geo 3 2 OUTPUT beam_3d_quad.msh) add_mesh(test_implicit_beam_3d_1 beam_3d.geo 3 1 OUTPUT beam_3d_lin.msh) register_test(test_solid_mechanics_model_implicit_dynamic_2d SOURCES test_solid_mechanics_model_implicit_dynamic_2d.cc DEPENDENCIES test_implicit_beam_2d_1 test_implicit_beam_2d_2 test_implicit_beam_3d_2 test_implicit_beam_3d_1 FILES_TO_COPY material_implicit_dynamic.dat DIRECTORIES_TO_CREATE paraview ) #=============================================================================== add_mesh(test_pbc_parallel_mesh square_structured.geo 2 1 OUTPUT square_structured.msh) register_test(test_solid_mechanics_model_bar_traction2d_structured_pbc SOURCES test_solid_mechanics_model_bar_traction2d_structured_pbc.cc DEPENDENCIES test_bar_traction_2d_mesh_structured1 FILES_TO_COPY material.dat test_cst_energy.pl DIRECTORIES_TO_CREATE paraview ) #register_test(test_solid_mechanics_model_pbc_parallel # SOURCES test_solid_mechanics_model_pbc_parallel.cc # DEPENDENCIES test_pbc_parallel_mesh # FILES_TO_COPY material.dat # DIRECTORIES_TO_CREATE paraview # ) #=============================================================================== add_mesh(test_cube3d_mesh1 cube.geo 3 1 OUTPUT cube1.msh) add_mesh(test_cube3d_mesh2 cube.geo 3 2 OUTPUT cube2.msh) add_mesh(test_cube3d_mesh_structured cube_structured.geo 3 1 OUTPUT cube_structured.msh) register_test(test_solid_mechanics_model_cube3d SOURCES test_solid_mechanics_model_cube3d.cc DEPENDENCIES test_cube3d_mesh1 FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) register_test(test_solid_mechanics_model_cube3d_tetra10 SOURCES test_solid_mechanics_model_cube3d_tetra10.cc DEPENDENCIES test_cube3d_mesh2 FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) register_test(test_solid_mechanics_model_cube3d_pbc SOURCES test_solid_mechanics_model_cube3d_pbc.cc DEPENDENCIES test_cube3d_mesh_structured FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) #add_mesh(test_solid_mechanics_model_boundary_condition_mesh cube_physical_names.geo 3 1) register_test(test_solid_mechanics_model_boundary_condition SOURCES test_solid_mechanics_model_boundary_condition.cc FILES_TO_COPY material.dat cube_physical_names.msh ) + +#=============================================================================== +add_mesh(test_cube3d_two_mat_mesh cube_two_materials.geo 3 1) + +register_test(test_solid_mechanics_model_reassign_material + SOURCES test_solid_mechanics_model_reassign_material.cc + DEPENDENCIES test_cube3d_two_mat_mesh + FILES_TO_COPY two_materials.dat + ) + +#=============================================================================== +register_test(test_solid_mechanics_model_prestrain_material + SOURCES test_solid_mechanics_model_prestrain_material.cc + FILES_TO_COPY cube_3d_tet_4.msh; material_elastic_plane_strain.dat + ) \ No newline at end of file diff --git a/test/test_model/test_solid_mechanics_model/cube_3d_tet_4.msh b/test/test_model/test_solid_mechanics_model/cube_3d_tet_4.msh new file mode 100644 index 000000000..6e1cb31b2 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/cube_3d_tet_4.msh @@ -0,0 +1,258 @@ +$MeshFormat +2.1 0 8 +$EndMeshFormat +$Nodes +42 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0 0 1 +6 1 0 1 +7 0 1 1 +8 1 1 1 +9 0.4999999999990217 1 1 +10 1 1 0.5000000000013878 +11 0.5000000000013878 1 0 +12 0 1 0.4999999999990217 +13 0 0 0.4999999999990217 +14 0.4999999999990217 0 1 +15 1 0 0.5000000000013878 +16 0.5000000000013878 0 0 +17 0 0.5000000000013878 0 +18 0 0.5000000000013878 1 +19 1 0.5000000000013878 1 +20 1 0.5000000000013878 0 +21 0.6428571428567235 0.5000000000003965 1 +22 0.2785714285710512 0.3500000000004956 1 +23 0.2842857142853593 0.670000000000456 1 +24 0.6840722533465502 1 0.2675114780935685 +25 0.5886704076322169 1 0.6013773593859361 +26 0.354548532196031 1 0.2737777674957053 +27 0.5000000000000585 0 0.6428571428572305 +28 0.3500000000004281 0 0.2785714285713892 +29 0.670000000000375 0 0.2842857142860015 +30 0.3571428571434519 0.5000000000003965 0 +31 0.7214285714291068 0.3500000000004956 0 +32 0.7157142857147892 0.670000000000456 0 +33 0 0.6428571428577375 0.4999999999997205 +34 0 0.2785714285719638 0.6499999999996506 +35 0 0.2842857142862178 0.3299999999996786 +36 1 0.3571428571434519 0.5000000000003965 +37 1 0.7214285714291068 0.3500000000004956 +38 1 0.7157142857147892 0.670000000000456 +39 0.4662426486062489 0.682365878613488 0.6247385391772754 +40 0.3928508996147502 0.6783805460075457 0.2231924003785042 +41 0.4834868113081788 0.3674115946766319 0.2900665800443959 +42 0.4329296945943543 0.3181780964219008 0.5652619629649539 +$EndNodes +$Elements +207 +1 15 3 0 1 0 1 +2 15 3 0 2 0 2 +3 15 3 0 3 0 3 +4 15 3 0 4 0 4 +5 15 3 0 5 0 5 +6 15 3 0 6 0 6 +7 15 3 0 7 0 7 +8 15 3 0 8 0 8 +9 1 3 0 1 0 7 9 +10 1 3 0 1 0 9 8 +11 1 3 0 2 0 8 10 +12 1 3 0 2 0 10 4 +13 1 3 0 3 0 4 11 +14 1 3 0 3 0 11 3 +15 1 3 0 4 0 3 12 +16 1 3 0 4 0 12 7 +17 1 3 0 5 0 1 13 +18 1 3 0 5 0 13 5 +19 1 3 0 6 0 5 14 +20 1 3 0 6 0 14 6 +21 1 3 0 7 0 6 15 +22 1 3 0 7 0 15 2 +23 1 3 0 8 0 2 16 +24 1 3 0 8 0 16 1 +25 1 3 0 9 0 3 17 +26 1 3 0 9 0 17 1 +27 1 3 0 10 0 7 18 +28 1 3 0 10 0 18 5 +29 1 3 0 11 0 8 19 +30 1 3 0 11 0 19 6 +31 1 3 0 12 0 4 20 +32 1 3 0 12 0 20 2 +33 2 3 0 14 0 21 9 8 +34 2 3 0 14 0 19 21 8 +35 2 3 0 14 0 14 21 6 +36 2 3 0 14 0 21 19 6 +37 2 3 0 14 0 22 14 5 +38 2 3 0 14 0 18 22 5 +39 2 3 0 14 0 14 22 21 +40 2 3 0 14 0 9 23 7 +41 2 3 0 14 0 23 18 7 +42 2 3 0 14 0 23 9 21 +43 2 3 0 14 0 23 22 18 +44 2 3 0 14 0 21 22 23 +45 2 3 0 16 0 24 4 11 +46 2 3 0 16 0 10 4 24 +47 2 3 0 16 0 25 8 10 +48 2 3 0 16 0 9 8 25 +49 2 3 0 16 0 25 10 24 +50 2 3 0 16 0 11 3 26 +51 2 3 0 16 0 26 3 12 +52 2 3 0 16 0 24 11 26 +53 2 3 0 16 0 25 24 26 +54 2 3 0 16 0 12 7 9 +55 2 3 0 16 0 9 25 12 +56 2 3 0 16 0 25 26 12 +57 2 3 0 18 0 15 27 6 +58 2 3 0 18 0 27 14 6 +59 2 3 0 18 0 27 13 5 +60 2 3 0 18 0 14 27 5 +61 2 3 0 18 0 13 28 1 +62 2 3 0 18 0 28 16 1 +63 2 3 0 18 0 28 13 27 +64 2 3 0 18 0 29 15 2 +65 2 3 0 18 0 16 29 2 +66 2 3 0 18 0 15 29 27 +67 2 3 0 18 0 28 29 16 +68 2 3 0 18 0 28 27 29 +69 2 3 0 20 0 11 3 30 +70 2 3 0 20 0 30 3 17 +71 2 3 0 20 0 30 1 16 +72 2 3 0 20 0 17 1 30 +73 2 3 0 20 0 16 2 31 +74 2 3 0 20 0 31 2 20 +75 2 3 0 20 0 31 30 16 +76 2 3 0 20 0 32 4 11 +77 2 3 0 20 0 20 4 32 +78 2 3 0 20 0 11 30 32 +79 2 3 0 20 0 31 20 32 +80 2 3 0 20 0 31 32 30 +81 2 3 0 22 0 33 17 3 +82 2 3 0 22 0 12 33 3 +83 2 3 0 22 0 33 12 7 +84 2 3 0 22 0 18 33 7 +85 2 3 0 22 0 34 18 5 +86 2 3 0 22 0 18 34 33 +87 2 3 0 22 0 13 34 5 +88 2 3 0 22 0 35 13 1 +89 2 3 0 22 0 35 17 33 +90 2 3 0 22 0 17 35 1 +91 2 3 0 22 0 35 34 13 +92 2 3 0 22 0 33 34 35 +93 2 3 0 24 0 36 19 6 +94 2 3 0 24 0 15 36 6 +95 2 3 0 24 0 36 15 2 +96 2 3 0 24 0 20 36 2 +97 2 3 0 24 0 37 20 4 +98 2 3 0 24 0 20 37 36 +99 2 3 0 24 0 10 37 4 +100 2 3 0 24 0 38 10 8 +101 2 3 0 24 0 38 19 36 +102 2 3 0 24 0 19 38 8 +103 2 3 0 24 0 38 37 10 +104 2 3 0 24 0 36 37 38 +105 4 3 0 26 0 38 25 37 39 +106 4 3 0 26 0 12 3 26 33 +107 4 3 0 26 0 4 20 37 32 +108 4 3 0 26 0 31 40 41 37 +109 4 3 0 26 0 28 35 41 1 +110 4 3 0 26 0 31 20 32 37 +111 4 3 0 26 0 40 3 26 11 +112 4 3 0 26 0 12 33 26 39 +113 4 3 0 26 0 9 12 25 39 +114 4 3 0 26 0 7 18 33 23 +115 4 3 0 26 0 42 33 34 39 +116 4 3 0 26 0 22 42 34 39 +117 4 3 0 26 0 42 29 36 27 +118 4 3 0 26 0 40 3 17 33 +119 4 3 0 26 0 21 36 6 27 +120 4 3 0 26 0 18 33 23 34 +121 4 3 0 26 0 22 18 23 34 +122 4 3 0 26 0 36 38 37 39 +123 4 3 0 26 0 22 21 39 23 +124 4 3 0 26 0 13 35 28 1 +125 4 3 0 26 0 5 22 34 18 +126 4 3 0 26 0 25 12 26 39 +127 4 3 0 26 0 30 3 17 40 +128 4 3 0 26 0 35 42 41 33 +129 4 3 0 26 0 22 34 23 39 +130 4 3 0 26 0 16 28 41 1 +131 4 3 0 26 0 21 38 9 8 +132 4 3 0 26 0 31 30 40 32 +133 4 3 0 26 0 36 29 2 15 +134 4 3 0 26 0 31 29 2 36 +135 4 3 0 26 0 35 40 17 33 +136 4 3 0 26 0 42 29 28 41 +137 4 3 0 26 0 36 31 41 37 +138 4 3 0 26 0 21 9 39 23 +139 4 3 0 26 0 13 5 27 34 +140 4 3 0 26 0 35 42 28 41 +141 4 3 0 26 0 31 30 41 40 +142 4 3 0 26 0 16 29 41 28 +143 4 3 0 26 0 24 10 37 4 +144 4 3 0 26 0 30 3 40 11 +145 4 3 0 26 0 42 36 41 39 +146 4 3 0 26 0 24 4 37 32 +147 4 3 0 26 0 40 35 41 33 +148 4 3 0 26 0 24 40 32 37 +149 4 3 0 26 0 16 31 41 29 +150 4 3 0 26 0 36 29 15 27 +151 4 3 0 26 0 21 38 19 36 +152 4 3 0 26 0 38 21 19 8 +153 4 3 0 26 0 12 7 39 9 +154 4 3 0 26 0 40 24 32 11 +155 4 3 0 26 0 24 40 26 11 +156 4 3 0 26 0 21 36 19 6 +157 4 3 0 26 0 4 24 11 32 +158 4 3 0 26 0 21 22 39 42 +159 4 3 0 26 0 20 31 2 36 +160 4 3 0 26 0 16 31 29 2 +161 4 3 0 26 0 21 42 36 27 +162 4 3 0 26 0 38 10 37 25 +163 4 3 0 26 0 10 24 37 25 +164 4 3 0 26 0 14 21 6 27 +165 4 3 0 26 0 20 31 36 37 +166 4 3 0 26 0 42 29 27 28 +167 4 3 0 26 0 10 38 8 25 +168 4 3 0 26 0 42 29 41 36 +169 4 3 0 26 0 42 35 34 33 +170 4 3 0 26 0 15 36 27 6 +171 4 3 0 26 0 34 33 23 39 +172 4 3 0 26 0 29 31 41 36 +173 4 3 0 26 0 31 16 41 30 +174 4 3 0 26 0 38 9 8 25 +175 4 3 0 26 0 30 40 32 11 +176 4 3 0 26 0 9 7 39 23 +177 4 3 0 26 0 30 16 41 1 +178 4 3 0 26 0 33 40 26 39 +179 4 3 0 26 0 33 42 41 39 +180 4 3 0 26 0 3 40 26 33 +181 4 3 0 26 0 40 33 41 39 +182 4 3 0 26 0 40 31 32 37 +183 4 3 0 26 0 27 34 28 13 +184 4 3 0 26 0 28 34 27 42 +185 4 3 0 26 0 34 35 28 13 +186 4 3 0 26 0 28 35 34 42 +187 4 3 0 26 0 30 1 35 17 +188 4 3 0 26 0 35 1 30 41 +189 4 3 0 26 0 30 35 40 17 +190 4 3 0 26 0 40 35 30 41 +191 4 3 0 26 0 41 37 39 36 +192 4 3 0 26 0 39 37 41 40 +193 4 3 0 26 0 36 39 21 38 +194 4 3 0 26 0 21 39 36 42 +195 4 3 0 26 0 26 39 24 25 +196 4 3 0 26 0 24 39 26 40 +197 4 3 0 26 0 39 37 24 25 +198 4 3 0 26 0 24 37 39 40 +199 4 3 0 26 0 7 39 33 12 +200 4 3 0 26 0 33 39 7 23 +201 4 3 0 26 0 27 21 22 14 +202 4 3 0 26 0 27 22 21 42 +203 4 3 0 26 0 34 22 27 42 +204 4 3 0 26 0 9 39 38 21 +205 4 3 0 26 0 38 39 9 25 +206 4 3 0 26 0 27 5 22 34 +207 4 3 0 26 0 27 22 5 14 +$EndElements diff --git a/test/test_model/test_solid_mechanics_model/cube_two_materials.geo b/test/test_model/test_solid_mechanics_model/cube_two_materials.geo new file mode 100644 index 000000000..2b9c23319 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/cube_two_materials.geo @@ -0,0 +1,14 @@ +L = 1; +s = 0.4; + +Point(1) = {-L/2., -L/2., 0, s}; + +line[] = Extrude{L,0,0}{ Point{1}; }; +surface[] = Extrude{0,L,0}{ Line{line[1]}; }; +volume[] = Extrude{0,0,L}{ Surface{surface[1]}; }; + +Physical Volume(1) = {volume[1]}; + +Transfinite Surface "*"; +Recombine Surface "*"; +Transfinite Volume "*"; diff --git a/test/test_model/test_solid_mechanics_model/material_elastic_plane_strain.dat b/test/test_model/test_solid_mechanics_model/material_elastic_plane_strain.dat new file mode 100644 index 000000000..078a85816 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/material_elastic_plane_strain.dat @@ -0,0 +1,8 @@ + +material elastic [ + name = steel + rho = 7800 # density + E = 2.1e11 # young's modulus + nu = 0.3 # poisson's ratio + Plane_Stress = 0 # plane strain +] diff --git a/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt index 82196d61d..84b468a37 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/test_materials/CMakeLists.txt @@ -1,61 +1,65 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date Fri Nov 26 00:17:56 2010 # # @brief configuration for materials tests # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_mesh(test_local_material_barre_trou_mesh barre_trou.geo 2 2) register_test(test_local_material SOURCES test_local_material.cc local_material_damage.cc EXTRA_FILES local_material_damage.hh local_material_damage_inline_impl.cc DEPENDENCIES test_local_material_barre_trou_mesh FILES_TO_COPY material.dat DIRECTORIES_TO_CREATE paraview ) add_mesh(test_material_thermal_mesh square.geo 2 1) register_test(test_material_thermal SOURCES test_material_thermal.cc DEPENDENCIES test_material_thermal_mesh FILES_TO_COPY material_thermal.dat ) #=============================================================================== add_mesh(test_interpolate_stress_mesh interpolation.geo 3 2) register_test(test_interpolate_stress test_interpolate_stress.cc FILES_TO_COPY material.dat DEPENDENCIES test_interpolate_stress_mesh DIRECTORIES_TO_CREATE paraview) #=============================================================================== + add_akantu_test(test_material_plasticity "test material plasticity" PACKAGE extra_materials) +add_akantu_test(test_material_damage "test material damage" + PACKAGE extra_materials) + add_akantu_test(test_material_viscoelastic "test the visco elastic materials" PACKAGE extra_materials) add_akantu_test(test_material_non_local "test the non-local materials" PACKAGE damage_non_local) diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/CMakeLists.txt new file mode 100644 index 000000000..b708e357c --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/CMakeLists.txt @@ -0,0 +1,37 @@ +#=============================================================================== +# @file CMakeLists.txt +# +# @author Nicolas Richart +# @author Aurelia Cuba Ramos +# +# @date Mon Jan 13 16:05:48 2014 +# +# @brief configuration for the material_damage tests +# +# @section LICENSE +# +# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . +# +# @section DESCRIPTION +# +#=============================================================================== + +register_test(test_material_damage + SOURCES test_material_damage.cc + FILES_TO_COPY cube_3d_tet_4.msh; material_damage.dat + PACKAGE extra_materials + ) diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/cube_3d_tet_4.msh b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/cube_3d_tet_4.msh new file mode 100644 index 000000000..6e1cb31b2 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/cube_3d_tet_4.msh @@ -0,0 +1,258 @@ +$MeshFormat +2.1 0 8 +$EndMeshFormat +$Nodes +42 +1 0 0 0 +2 1 0 0 +3 0 1 0 +4 1 1 0 +5 0 0 1 +6 1 0 1 +7 0 1 1 +8 1 1 1 +9 0.4999999999990217 1 1 +10 1 1 0.5000000000013878 +11 0.5000000000013878 1 0 +12 0 1 0.4999999999990217 +13 0 0 0.4999999999990217 +14 0.4999999999990217 0 1 +15 1 0 0.5000000000013878 +16 0.5000000000013878 0 0 +17 0 0.5000000000013878 0 +18 0 0.5000000000013878 1 +19 1 0.5000000000013878 1 +20 1 0.5000000000013878 0 +21 0.6428571428567235 0.5000000000003965 1 +22 0.2785714285710512 0.3500000000004956 1 +23 0.2842857142853593 0.670000000000456 1 +24 0.6840722533465502 1 0.2675114780935685 +25 0.5886704076322169 1 0.6013773593859361 +26 0.354548532196031 1 0.2737777674957053 +27 0.5000000000000585 0 0.6428571428572305 +28 0.3500000000004281 0 0.2785714285713892 +29 0.670000000000375 0 0.2842857142860015 +30 0.3571428571434519 0.5000000000003965 0 +31 0.7214285714291068 0.3500000000004956 0 +32 0.7157142857147892 0.670000000000456 0 +33 0 0.6428571428577375 0.4999999999997205 +34 0 0.2785714285719638 0.6499999999996506 +35 0 0.2842857142862178 0.3299999999996786 +36 1 0.3571428571434519 0.5000000000003965 +37 1 0.7214285714291068 0.3500000000004956 +38 1 0.7157142857147892 0.670000000000456 +39 0.4662426486062489 0.682365878613488 0.6247385391772754 +40 0.3928508996147502 0.6783805460075457 0.2231924003785042 +41 0.4834868113081788 0.3674115946766319 0.2900665800443959 +42 0.4329296945943543 0.3181780964219008 0.5652619629649539 +$EndNodes +$Elements +207 +1 15 3 0 1 0 1 +2 15 3 0 2 0 2 +3 15 3 0 3 0 3 +4 15 3 0 4 0 4 +5 15 3 0 5 0 5 +6 15 3 0 6 0 6 +7 15 3 0 7 0 7 +8 15 3 0 8 0 8 +9 1 3 0 1 0 7 9 +10 1 3 0 1 0 9 8 +11 1 3 0 2 0 8 10 +12 1 3 0 2 0 10 4 +13 1 3 0 3 0 4 11 +14 1 3 0 3 0 11 3 +15 1 3 0 4 0 3 12 +16 1 3 0 4 0 12 7 +17 1 3 0 5 0 1 13 +18 1 3 0 5 0 13 5 +19 1 3 0 6 0 5 14 +20 1 3 0 6 0 14 6 +21 1 3 0 7 0 6 15 +22 1 3 0 7 0 15 2 +23 1 3 0 8 0 2 16 +24 1 3 0 8 0 16 1 +25 1 3 0 9 0 3 17 +26 1 3 0 9 0 17 1 +27 1 3 0 10 0 7 18 +28 1 3 0 10 0 18 5 +29 1 3 0 11 0 8 19 +30 1 3 0 11 0 19 6 +31 1 3 0 12 0 4 20 +32 1 3 0 12 0 20 2 +33 2 3 0 14 0 21 9 8 +34 2 3 0 14 0 19 21 8 +35 2 3 0 14 0 14 21 6 +36 2 3 0 14 0 21 19 6 +37 2 3 0 14 0 22 14 5 +38 2 3 0 14 0 18 22 5 +39 2 3 0 14 0 14 22 21 +40 2 3 0 14 0 9 23 7 +41 2 3 0 14 0 23 18 7 +42 2 3 0 14 0 23 9 21 +43 2 3 0 14 0 23 22 18 +44 2 3 0 14 0 21 22 23 +45 2 3 0 16 0 24 4 11 +46 2 3 0 16 0 10 4 24 +47 2 3 0 16 0 25 8 10 +48 2 3 0 16 0 9 8 25 +49 2 3 0 16 0 25 10 24 +50 2 3 0 16 0 11 3 26 +51 2 3 0 16 0 26 3 12 +52 2 3 0 16 0 24 11 26 +53 2 3 0 16 0 25 24 26 +54 2 3 0 16 0 12 7 9 +55 2 3 0 16 0 9 25 12 +56 2 3 0 16 0 25 26 12 +57 2 3 0 18 0 15 27 6 +58 2 3 0 18 0 27 14 6 +59 2 3 0 18 0 27 13 5 +60 2 3 0 18 0 14 27 5 +61 2 3 0 18 0 13 28 1 +62 2 3 0 18 0 28 16 1 +63 2 3 0 18 0 28 13 27 +64 2 3 0 18 0 29 15 2 +65 2 3 0 18 0 16 29 2 +66 2 3 0 18 0 15 29 27 +67 2 3 0 18 0 28 29 16 +68 2 3 0 18 0 28 27 29 +69 2 3 0 20 0 11 3 30 +70 2 3 0 20 0 30 3 17 +71 2 3 0 20 0 30 1 16 +72 2 3 0 20 0 17 1 30 +73 2 3 0 20 0 16 2 31 +74 2 3 0 20 0 31 2 20 +75 2 3 0 20 0 31 30 16 +76 2 3 0 20 0 32 4 11 +77 2 3 0 20 0 20 4 32 +78 2 3 0 20 0 11 30 32 +79 2 3 0 20 0 31 20 32 +80 2 3 0 20 0 31 32 30 +81 2 3 0 22 0 33 17 3 +82 2 3 0 22 0 12 33 3 +83 2 3 0 22 0 33 12 7 +84 2 3 0 22 0 18 33 7 +85 2 3 0 22 0 34 18 5 +86 2 3 0 22 0 18 34 33 +87 2 3 0 22 0 13 34 5 +88 2 3 0 22 0 35 13 1 +89 2 3 0 22 0 35 17 33 +90 2 3 0 22 0 17 35 1 +91 2 3 0 22 0 35 34 13 +92 2 3 0 22 0 33 34 35 +93 2 3 0 24 0 36 19 6 +94 2 3 0 24 0 15 36 6 +95 2 3 0 24 0 36 15 2 +96 2 3 0 24 0 20 36 2 +97 2 3 0 24 0 37 20 4 +98 2 3 0 24 0 20 37 36 +99 2 3 0 24 0 10 37 4 +100 2 3 0 24 0 38 10 8 +101 2 3 0 24 0 38 19 36 +102 2 3 0 24 0 19 38 8 +103 2 3 0 24 0 38 37 10 +104 2 3 0 24 0 36 37 38 +105 4 3 0 26 0 38 25 37 39 +106 4 3 0 26 0 12 3 26 33 +107 4 3 0 26 0 4 20 37 32 +108 4 3 0 26 0 31 40 41 37 +109 4 3 0 26 0 28 35 41 1 +110 4 3 0 26 0 31 20 32 37 +111 4 3 0 26 0 40 3 26 11 +112 4 3 0 26 0 12 33 26 39 +113 4 3 0 26 0 9 12 25 39 +114 4 3 0 26 0 7 18 33 23 +115 4 3 0 26 0 42 33 34 39 +116 4 3 0 26 0 22 42 34 39 +117 4 3 0 26 0 42 29 36 27 +118 4 3 0 26 0 40 3 17 33 +119 4 3 0 26 0 21 36 6 27 +120 4 3 0 26 0 18 33 23 34 +121 4 3 0 26 0 22 18 23 34 +122 4 3 0 26 0 36 38 37 39 +123 4 3 0 26 0 22 21 39 23 +124 4 3 0 26 0 13 35 28 1 +125 4 3 0 26 0 5 22 34 18 +126 4 3 0 26 0 25 12 26 39 +127 4 3 0 26 0 30 3 17 40 +128 4 3 0 26 0 35 42 41 33 +129 4 3 0 26 0 22 34 23 39 +130 4 3 0 26 0 16 28 41 1 +131 4 3 0 26 0 21 38 9 8 +132 4 3 0 26 0 31 30 40 32 +133 4 3 0 26 0 36 29 2 15 +134 4 3 0 26 0 31 29 2 36 +135 4 3 0 26 0 35 40 17 33 +136 4 3 0 26 0 42 29 28 41 +137 4 3 0 26 0 36 31 41 37 +138 4 3 0 26 0 21 9 39 23 +139 4 3 0 26 0 13 5 27 34 +140 4 3 0 26 0 35 42 28 41 +141 4 3 0 26 0 31 30 41 40 +142 4 3 0 26 0 16 29 41 28 +143 4 3 0 26 0 24 10 37 4 +144 4 3 0 26 0 30 3 40 11 +145 4 3 0 26 0 42 36 41 39 +146 4 3 0 26 0 24 4 37 32 +147 4 3 0 26 0 40 35 41 33 +148 4 3 0 26 0 24 40 32 37 +149 4 3 0 26 0 16 31 41 29 +150 4 3 0 26 0 36 29 15 27 +151 4 3 0 26 0 21 38 19 36 +152 4 3 0 26 0 38 21 19 8 +153 4 3 0 26 0 12 7 39 9 +154 4 3 0 26 0 40 24 32 11 +155 4 3 0 26 0 24 40 26 11 +156 4 3 0 26 0 21 36 19 6 +157 4 3 0 26 0 4 24 11 32 +158 4 3 0 26 0 21 22 39 42 +159 4 3 0 26 0 20 31 2 36 +160 4 3 0 26 0 16 31 29 2 +161 4 3 0 26 0 21 42 36 27 +162 4 3 0 26 0 38 10 37 25 +163 4 3 0 26 0 10 24 37 25 +164 4 3 0 26 0 14 21 6 27 +165 4 3 0 26 0 20 31 36 37 +166 4 3 0 26 0 42 29 27 28 +167 4 3 0 26 0 10 38 8 25 +168 4 3 0 26 0 42 29 41 36 +169 4 3 0 26 0 42 35 34 33 +170 4 3 0 26 0 15 36 27 6 +171 4 3 0 26 0 34 33 23 39 +172 4 3 0 26 0 29 31 41 36 +173 4 3 0 26 0 31 16 41 30 +174 4 3 0 26 0 38 9 8 25 +175 4 3 0 26 0 30 40 32 11 +176 4 3 0 26 0 9 7 39 23 +177 4 3 0 26 0 30 16 41 1 +178 4 3 0 26 0 33 40 26 39 +179 4 3 0 26 0 33 42 41 39 +180 4 3 0 26 0 3 40 26 33 +181 4 3 0 26 0 40 33 41 39 +182 4 3 0 26 0 40 31 32 37 +183 4 3 0 26 0 27 34 28 13 +184 4 3 0 26 0 28 34 27 42 +185 4 3 0 26 0 34 35 28 13 +186 4 3 0 26 0 28 35 34 42 +187 4 3 0 26 0 30 1 35 17 +188 4 3 0 26 0 35 1 30 41 +189 4 3 0 26 0 30 35 40 17 +190 4 3 0 26 0 40 35 30 41 +191 4 3 0 26 0 41 37 39 36 +192 4 3 0 26 0 39 37 41 40 +193 4 3 0 26 0 36 39 21 38 +194 4 3 0 26 0 21 39 36 42 +195 4 3 0 26 0 26 39 24 25 +196 4 3 0 26 0 24 39 26 40 +197 4 3 0 26 0 39 37 24 25 +198 4 3 0 26 0 24 37 39 40 +199 4 3 0 26 0 7 39 33 12 +200 4 3 0 26 0 33 39 7 23 +201 4 3 0 26 0 27 21 22 14 +202 4 3 0 26 0 27 22 21 42 +203 4 3 0 26 0 34 22 27 42 +204 4 3 0 26 0 9 39 38 21 +205 4 3 0 26 0 38 39 9 25 +206 4 3 0 26 0 27 5 22 34 +207 4 3 0 26 0 27 22 5 14 +$EndElements diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/material_damage.dat b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/material_damage.dat new file mode 100644 index 000000000..94313d906 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/material_damage.dat @@ -0,0 +1,8 @@ + +material damage [ + name = steel + rho = 7800 # density + E = 2.1e11 # young's modulus + nu = 0.3 # poisson's ratio + Plane_Stress = 0 # plane strain +] diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/test_material_damage.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/test_material_damage.cc new file mode 100644 index 000000000..87423c44d --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage/test_material_damage.cc @@ -0,0 +1,233 @@ +/** + * @file test_material_damage.cc + * + * @author Nicolas Richart + * @author Aurelia Cuba Ramos + * + * @date Mon Jan 13 16:05:48 2014 + * + * @brief test for material damage class + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ + +#include "aka_common.hh" +#include "mesh.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model.hh" +#include "material.hh" +#include "element_class.hh" + +using namespace akantu; + +Real alpha [3][4] = { { 0.01, 0.02, 0.03, 0.04 }, + { 0.05, 0.06, 0.07, 0.08 }, + { 0.09, 0.10, 0.11, 0.12 } }; + +/* -------------------------------------------------------------------------- */ +template +static Matrix prescribed_strain() { + UInt spatial_dimension = ElementClass::getSpatialDimension(); + Matrix strain(spatial_dimension, spatial_dimension); + + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { + strain(i, j) = alpha[i][j + 1]; + } + } + return strain; +} + +template +static Matrix prescribed_stress(Real prescribed_dam) { + UInt spatial_dimension = ElementClass::getSpatialDimension(); + Matrix stress(spatial_dimension, spatial_dimension); + + //plane strain in 2d + Matrix strain(spatial_dimension, spatial_dimension); + Matrix pstrain; pstrain = prescribed_strain(); + Real nu = 0.3; + Real E = 2.1e11; + Real trace = 0; + + /// symetric part of the strain tensor + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); + + + for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); + + Real lambda = nu * E / ((1 + nu) * (1 - 2*nu)); + Real mu = E / (2 * (1 + nu)); + + if(!is_plane_strain) { + std::cout << "toto" << std::endl; + lambda = nu * E / (1 - nu*nu); + } + + if(spatial_dimension == 1) { + stress(0, 0) = (1-prescribed_dam) * E * strain(0, 0); + } else { + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) { + stress(i, j) = (i == j)*(1-prescribed_dam)*lambda*trace + 2*mu*(1-prescribed_dam)*strain(i, j); + } + } + + return stress; +} + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) +{ + initialize(argc, argv); + + Real prescribed_dam = 0.1; + + UInt dim = 3; + const ElementType element_type = _tetrahedron_4; + const bool plane_strain = true; + + /// load mesh + Mesh my_mesh(dim); + MeshIOMSH mesh_io; + + std::stringstream filename; filename << "cube_3d_tet_4.msh"; + mesh_io.read(filename.str(), my_mesh); + + // MeshUtils::purifyMesh(my_mesh); + + UInt nb_nodes = my_mesh.getNbNodes(); + + /// declaration of model + SolidMechanicsModel my_model(my_mesh); + /// model initialization + my_model.initFull("material_damage.dat", SolidMechanicsModelOptions(_static)); + + const Array & coordinates = my_mesh.getNodes(); + Array & displacement = my_model.getDisplacement(); + Array & boundary = my_model.getBoundary(); + MeshUtils::buildFacets(my_mesh); + + my_mesh.createBoundaryGroupFromGeometry(); + + // Loop over (Sub)Boundar(ies) + for(GroupManager::const_element_group_iterator it(my_mesh.element_group_begin()); + it != my_mesh.element_group_end(); ++it) { + for(ElementGroup::const_node_iterator nodes_it(it->second->node_begin()); + nodes_it!= it->second->node_end(); ++nodes_it) { + UInt n(*nodes_it); + std::cout << "Node " << *nodes_it << std::endl; + for (UInt i = 0; i < dim; ++i) { + displacement(n, i) = alpha[i][0]; + for (UInt j = 0; j < dim; ++j) { + displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); + } + boundary(n, i) = true; + } + } + } + + /* ------------------------------------------------------------------------ */ + /* Set damage in each gauss point */ + /* ------------------------------------------------------------------------ */ + UInt nb_quadrature_points = my_model.getFEM().getNbQuadraturePoints(element_type); + + Array & dam_vect = const_cast &>(my_model.getMaterial(0).getInternal("damage")(element_type)); + Array::iterator dam_it = dam_vect.begin(); + + UInt nb_element = my_mesh.getNbElement(element_type); + + for (UInt el = 0; el < nb_element; ++el) { + for (UInt q= 0; q < nb_quadrature_points; ++q) { + *dam_it = prescribed_dam; + + ++dam_it; + } + } + /* ------------------------------------------------------------------------ */ + /* Static solve */ + /* ------------------------------------------------------------------------ */ + my_model.solveStep<_scm_newton_raphson_tangent_modified, _scc_residual>(2e-4, 2); + + /* ------------------------------------------------------------------------ */ + /* Checks */ + /* ------------------------------------------------------------------------ */ + + + Array & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); + + Array::iterator< Matrix > stress_it = stress_vect.begin(dim, dim); + + Matrix presc_stress; + presc_stress = prescribed_stress(prescribed_dam); + + Real stress_tolerance = 1e-13; + + for (UInt el = 0; el < nb_element; ++el) { + for (UInt q = 0; q < nb_quadrature_points; ++q) { + Matrix & stress = *stress_it; + + Matrix diff(dim, dim); + + diff = stress; + diff -= presc_stress; + Real stress_error = diff.norm() / stress.norm(); + + if(stress_error > stress_tolerance) { + std::cerr << "stress error: " << stress_error << " > " << stress_tolerance << std::endl; + std::cerr << "stress: " << stress << std::endl + << "prescribed stress: " << presc_stress << std::endl; + return EXIT_FAILURE; + } else { + std::cerr << "stress error: " << stress_error << " < " << stress_tolerance << std::endl; + } + + ++stress_it; + } + } + + + for (UInt n = 0; n < nb_nodes; ++n) { + for (UInt i = 0; i < dim; ++i) { + Real disp = alpha[i][0]; + for (UInt j = 0; j < dim; ++j) { + disp += alpha[i][j + 1] * coordinates(n, j); + } + + if(!(std::abs(displacement(n,i) - disp) < 2e-15)) { + std::cerr << "displacement(" << n << ", " << i <<")=" << displacement(n,i) << " should be equal to " << disp << std::endl; + return EXIT_FAILURE; + } + } + } + + finalize(); + + return EXIT_SUCCESS; +} + + diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_prestrain_material.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_prestrain_material.cc new file mode 100644 index 000000000..cdb907bd9 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_prestrain_material.cc @@ -0,0 +1,216 @@ +/** + * @file test_solid_mechanics_model_prestrain_material.cc + * + * @author Nicolas Richart + * @author Aurelia Cuba Ramos + * + * @date Mon Jan 13 16:05:48 2014 + * + * @brief test the internal field prestrain + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ + +#include "aka_common.hh" +#include "mesh.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model.hh" +#include "material.hh" +#include "element_class.hh" + +using namespace akantu; + +Real alpha [3][4] = { { 0.01, 0.02, 0.03, 0.04 }, + { 0.05, 0.06, 0.07, 0.08 }, + { 0.09, 0.10, 0.11, 0.12 } }; + +/* -------------------------------------------------------------------------- */ +template +static Matrix prescribed_strain() { + UInt spatial_dimension = ElementClass::getSpatialDimension(); + Matrix strain(spatial_dimension, spatial_dimension); + + for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { + strain(i, j) = alpha[i][j + 1]; + } + } + return strain; +} + +template +static Matrix prescribed_stress(Matrix prescribed_prestrain) { + UInt spatial_dimension = ElementClass::getSpatialDimension(); + Matrix stress(spatial_dimension, spatial_dimension); + + //plane strain in 2d + Matrix strain(spatial_dimension, spatial_dimension); + Matrix pstrain; + pstrain = prescribed_strain() + prescribed_prestrain; + Real nu = 0.3; + Real E = 2.1e11; + Real trace = 0; + + /// symetric part of the strain tensor + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + strain(i,j) = 0.5 * (pstrain(i, j) + pstrain(j, i)); + + + for (UInt i = 0; i < spatial_dimension; ++i) trace += strain(i,i); + + Real lambda = nu * E / ((1 + nu) * (1 - 2*nu)); + Real mu = E / (2 * (1 + nu)); + + if(!is_plane_strain) { + std::cout << "toto" << std::endl; + lambda = nu * E / (1 - nu*nu); + } + + if(spatial_dimension == 1) { + stress(0, 0) = E * strain(0, 0); + } else { + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) { + stress(i, j) = (i == j)*lambda*trace + 2*mu*strain(i, j); + } + } + + return stress; +} + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) +{ + initialize(argc, argv); + + UInt dim = 3; + const ElementType element_type = _tetrahedron_4; + const bool plane_strain = true; + Matrix prescribed_prestrain(dim, dim); + prescribed_prestrain.clear(); + for (UInt i = 0; i < dim; ++i) { + for (UInt j = 0; j < dim; ++j) + prescribed_prestrain(i,j) += 0.1; + } + + + /// load mesh + Mesh my_mesh(dim); + MeshIOMSH mesh_io; + + std::stringstream filename; filename << "cube_3d_tet_4.msh"; + mesh_io.read(filename.str(), my_mesh); + + /// declaration of model + SolidMechanicsModel my_model(my_mesh); + /// model initialization + my_model.initFull("material_elastic_plane_strain.dat", SolidMechanicsModelOptions(_static)); + + const Array & coordinates = my_mesh.getNodes(); + Array & displacement = my_model.getDisplacement(); + Array & boundary = my_model.getBoundary(); + MeshUtils::buildFacets(my_mesh); + + my_mesh.createBoundaryGroupFromGeometry(); + + // Loop over (Sub)Boundar(ies) + for(GroupManager::const_element_group_iterator it(my_mesh.element_group_begin()); + it != my_mesh.element_group_end(); ++it) { + for(ElementGroup::const_node_iterator nodes_it(it->second->node_begin()); + nodes_it!= it->second->node_end(); ++nodes_it) { + UInt n(*nodes_it); + std::cout << "Node " << *nodes_it << std::endl; + for (UInt i = 0; i < dim; ++i) { + displacement(n, i) = alpha[i][0]; + for (UInt j = 0; j < dim; ++j) { + displacement(n, i) += alpha[i][j + 1] * coordinates(n, j); + } + boundary(n, i) = true; + } + } + } + + /* ------------------------------------------------------------------------ */ + /* Apply prestrain in each element */ + /* ------------------------------------------------------------------------ */ + + + Array & prestrain_vect = const_cast &>(my_model.getMaterial(0).getInternal("pre_strain")(element_type)); + Array::iterator< Matrix > prestrain_it = prestrain_vect.begin(dim, dim); + Array::iterator< Matrix > prestrain_end = prestrain_vect.end(dim, dim); + + for (; prestrain_it != prestrain_end; ++prestrain_it) { + for (UInt i = 0; i < dim; ++i) + for (UInt j = 0; j < dim; ++j) + (*prestrain_it)(i,j) += prescribed_prestrain(i,j); + } + /* ------------------------------------------------------------------------ */ + /* Static solve */ + /* ------------------------------------------------------------------------ */ + my_model.solveStep<_scm_newton_raphson_tangent_modified, _scc_residual>(2e-4, 2); + + /* ------------------------------------------------------------------------ */ + /* Checks */ + /* ------------------------------------------------------------------------ */ + + + Array & stress_vect = const_cast &>(my_model.getMaterial(0).getStress(element_type)); + + Array::iterator< Matrix > stress_it = stress_vect.begin(dim, dim); + Array::iterator< Matrix > stress_end = stress_vect.end(dim, dim); + + Matrix presc_stress; + presc_stress = prescribed_stress(prescribed_prestrain); + + Real stress_tolerance = 1e-13; + + for (; stress_it != stress_end; ++stress_it) { + Matrix & stress = *stress_it; + Matrix diff(dim, dim); + + diff = stress; + diff -= presc_stress; + Real stress_error = diff.norm() / stress.norm(); + + if(stress_error > stress_tolerance) { + std::cerr << "stress error: " << stress_error << " > " << stress_tolerance << std::endl; + std::cerr << "stress: " << stress << std::endl + << "prescribed stress: " << presc_stress << std::endl; + return EXIT_FAILURE; + } else { + std::cerr << "stress error: " << stress_error << " < " << stress_tolerance << std::endl; + } + + + } + + + finalize(); + + return EXIT_SUCCESS; +} + + diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc new file mode 100644 index 000000000..bb4f868ab --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.cc @@ -0,0 +1,216 @@ +/** + * @file test_solid_mechanics_model_reassign_material.cc + * @author Aurelia Cuba Ramos + * @date Mon May 6 17:37:17 2013 + * + * @brief test the function reassign material + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "static_communicator.hh" +#include "solid_mechanics_model.hh" +#include "material.hh" +#include "aka_grid_dynamic.hh" +using namespace akantu; + + +class StraightInterfaceMaterialSelector : public MaterialSelector { +public: + StraightInterfaceMaterialSelector(SolidMechanicsModel & model, + const std::string & mat_1_material, + const std::string & mat_2_material, + bool & horizontal, + Real & pos_interface) : + model(model), + mat_1_material(mat_1_material), + mat_2_material(mat_2_material), + horizontal(horizontal), + pos_interface(pos_interface) { + Mesh & mesh = model.getMesh(); + UInt spatial_dimension = mesh.getSpatialDimension(); + + /// store barycenters of all elements + mesh.initByElementTypeArray(barycenters, spatial_dimension, spatial_dimension); + + for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { + GhostType ghost_type = *gt; + Element e; + e.ghost_type = ghost_type; + + Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); + for(; it != last_type; ++it) { + UInt nb_element = mesh.getNbElement(*it, ghost_type); + e.type = *it; + Array & barycenter = barycenters(*it, ghost_type); + barycenter.resize(nb_element); + + Array::iterator< Vector > bary_it = barycenter.begin(spatial_dimension); + for (UInt elem = 0; elem < nb_element; ++elem) { + e.element = elem; + mesh.getBarycenter(e, *bary_it); + ++bary_it; + } + } + } + } + + UInt operator()(const Element & elem) { + UInt spatial_dimension = model.getSpatialDimension(); + const Vector & bary = barycenters(elem.type, elem.ghost_type).begin(spatial_dimension)[elem.element]; + + /// check for a given element on which side of the material interface plane the bary center lies and assign corresponding material + if (bary(!horizontal) < pos_interface) { + return model.getMaterialIndex(mat_1_material);; + } + return model.getMaterialIndex(mat_2_material);; + + } + + bool isConditonVerified() { + + /// check if material has been (re)-assigned correctly + Mesh & mesh = model.getMesh(); + UInt spatial_dimension = mesh.getSpatialDimension(); + GhostType ghost_type = _not_ghost; + + Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); + for(; it != last_type; ++it) { + Array & el_idx_by_mat = model.getElementIndexByMaterial(*it, ghost_type); + UInt nb_element = mesh.getNbElement(*it, ghost_type); + Array::iterator > bary = barycenters(*it, ghost_type).begin(spatial_dimension); + for (UInt elem = 0; elem < nb_element; ++elem, ++bary) { + /// compare element_index_by material to material index that should be assigned due to the geometry of the interface + UInt mat_index; + if ((*bary)(!horizontal) < pos_interface) + mat_index = model.getMaterialIndex(mat_1_material); + else + mat_index = model.getMaterialIndex(mat_2_material); + + if (el_idx_by_mat(elem,0) != mat_index) + /// wrong material index, make test fail + return false; + + + } + } + return true; + } + + void moveInterface(Real & pos_new, bool & horizontal_new) { + /// update position and orientation of material interface plane + pos_interface = pos_new; + horizontal = horizontal_new; + model.reassignMaterial(); + } + +protected: + SolidMechanicsModel & model; + ByElementTypeReal barycenters; + std::string mat_1_material; + std::string mat_2_material; + bool horizontal; + Real pos_interface; +}; + + +/* -------------------------------------------------------------------------- */ +/* Main */ +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + + bool test_passed; + + debug::setDebugLevel(dblWarning); + initialize(argc, argv); + + /// specify position and orientation of material interface plane + bool horizontal = true; + Real pos_interface = 0.; + + UInt spatial_dimension = 3; + + akantu::StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); + akantu::Int psize = comm.getNbProc(); + akantu::Int prank = comm.whoAmI(); + + Mesh mesh(spatial_dimension); + + akantu::MeshPartition * partition = NULL; + if(prank == 0) { + /// creation mesh + mesh.read("cube_two_materials.msh"); + + partition = new akantu::MeshPartitionScotch(mesh, spatial_dimension); + partition->partitionate(psize); + } + + /// model creation + SolidMechanicsModel model(mesh); + model.initParallel(partition); + delete partition; + + /// assign the two different materials using the StraightInterfaceMaterialSelector + StraightInterfaceMaterialSelector * mat_selector; + mat_selector = new StraightInterfaceMaterialSelector(model, "mat_1", "mat_2", horizontal, pos_interface); + + model.setMaterialSelector(*mat_selector); + model.initFull("two_materials.dat", SolidMechanicsModelOptions(_static)); + + // model.setBaseName("test_reassign_material"); + // model.addDumpField("element_index_by_material"); + // model.addDumpField("partitions"); + // model.dump(); + + /// check if different materials have been assigned correctly + test_passed = mat_selector->isConditonVerified(); + if (!test_passed) { + AKANTU_DEBUG_ERROR("materials not correctly assigned"); + return EXIT_FAILURE; + } + + + + /// change orientation of material interface plane + horizontal = false; + mat_selector->moveInterface(pos_interface, horizontal); + + // model.dump(); + + /// test if material has been reassigned correctly + test_passed = mat_selector->isConditonVerified(); + if (!test_passed) { + AKANTU_DEBUG_ERROR("materials not correctly reassigned"); + return EXIT_FAILURE; + } + + finalize(); + + if(prank == 0) + std::cout << "OK: test passed!" << std::endl; + + return EXIT_SUCCESS; +} + +/* -------------------------------------------------------------------------- */ + diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.sh b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.sh new file mode 100755 index 000000000..7130b5154 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_reassign_material.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 8 ./test_solid_mechanics_model_reassign_material \ No newline at end of file diff --git a/test/test_model/test_solid_mechanics_model/two_materials.dat b/test/test_model/test_solid_mechanics_model/two_materials.dat new file mode 100644 index 000000000..2d15402fd --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/two_materials.dat @@ -0,0 +1,13 @@ +material elastic [ + name = mat_1 + rho = 2.2e-9 # density + E = 12e9 # young's modulus + nu = 0.3 # poisson's ratio +] + +material elastic [ + name = mat_2 + rho = 2.7e-9 # density + E = 60e9 # young's modulus + nu = 0.3 # poisson's ratio +]