diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index ec0b853e0..21a6c0a00 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,696 +1,780 @@ /** * @file material.cc * @author Nicolas Richart * @date Tue Jul 27 11:43:41 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(model.getMemoryID()), id(id), name(""), model(&model), stress("stress", id), strain("strain", id), element_filter("element_filter", id), // potential_energy_vector(false), potential_energy("potential_energy", id), - is_non_local(false), is_init(false) { + is_non_local(false), + interpolationInvCoordinates("interpolationInvCoordinates", id), + is_init(false) { AKANTU_DEBUG_IN(); rho = 0; AKANTU_DEBUG_ASSERT(this->model,"model has wrong type: cannot proceed"); spatial_dimension = this->model->getSpatialDimension(); /// allocate strain stress for local elements initInternalVector(strain, spatial_dimension * spatial_dimension); initInternalVector(stress, spatial_dimension * spatial_dimension); /// for each connectivity types allocate the element filer array of the material initInternalVector(element_filter, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool Material::setParam(const std::string & key, const std::string & value, __attribute__ ((unused)) const ID & id) { std::stringstream sstr(value); if(key == "name") name = std::string(value); else if(key == "rho") { sstr >> rho; } else return false; return true; } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const { AKANTU_DEBUG_IN(); model->getFEM().getMesh().initByElementTypeVector(vect, nb_component, spatial_dimension, false, element_kind); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::resizeInternalVector(ByElementTypeVector & by_el_type_vect, ElementKind element_kind) const { AKANTU_DEBUG_IN(); FEM * fem = & model->getFEM(); if (element_kind == _ek_cohesive) fem = & model->getFEM("CohesiveFEM"); const Mesh & mesh = fem->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, element_kind); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, element_kind); for(; it != end; ++it) { const Vector & elem_filter = element_filter(*it, gt); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, gt); UInt new_size = nb_element * nb_quadrature_points; Vector & vect = by_el_type_vect(*it, gt); UInt size = vect.getSize(); UInt nb_component = vect.getNbComponent(); vect.resize(new_size); memset(vect.storage() + size * nb_component, 0, (new_size - size) * nb_component * sizeof(T)); } } 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(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); computeStress(displacement, ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & 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 Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type); Vector & 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$ Vector * sigma_dphi_dx = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Real * shapesd = shapes_derivatives.storage(); Real * shapesd_val; UInt * elem_filter_val = elem_filter.storage(); Vector * shapesd_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd"); Real * shapesd_filtered_val = shapesd_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapesd_val = shapesd + elem_filter_val[el] * size_of_shapes_derivatives * nb_quadrature_points; memcpy(shapesd_filtered_val, shapesd_val, size_of_shapes_derivatives * nb_quadrature_points * sizeof(Real)); shapesd_filtered_val += size_of_shapes_derivatives * nb_quadrature_points; } Vector & stress_vect = stress(*it, ghost_type); // Vector::iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); // Vector::iterator sigma_end = stress_vect.end(spatial_dimension, spatial_dimension); // Vector::iterator nabla_B = shapesd_filtered->begin(nb_nodes_per_element, spatial_dimension); // Vector::iterator sigma_dphi_dx_it = sigma_dphi_dx->begin(nb_nodes_per_element, spatial_dimension); // for (; sigma != sigma_end; ++sigma, ++nabla_B, ++sigma_dphi_dx_it) { // sigma_dphi_dx_it->mul(*nabla_B, *sigma); // } Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, *shapesd_filtered, stress_vect, *sigma_dphi_dx); 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$ */ Vector * int_sigma_dphi_dx = new Vector(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().assembleVector(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, -1); delete int_sigma_dphi_dx; } 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::computeStress(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(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) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, 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::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Vector & displacement = model->getDisplacement(); resizeInternalVector(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) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & 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(Vector & current_position, 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: { assembleStiffnessMatrix<1>(current_position, *it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(current_position, *it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(current_position, *it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(type,ghost_type); Vector & elem_filter = element_filter(type, ghost_type); Vector & strain_vect = strain(type, ghost_type); UInt * elem_filter_val = elem_filter.storage(); UInt nb_element = elem_filter.getSize(); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); 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(current_position, strain_vect, dim, type, ghost_type, &elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Vector * tangent_stiffness_matrix = new Vector(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); computeTangentStiffness(type, *tangent_stiffness_matrix, ghost_type); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Vector * bt_d_b = new Vector(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); UInt size_of_b = tangent_size * bt_d_b_size; Real * B = new Real[size_of_b]; Real * Bt_D = new Real[size_of_b]; Real * Bt_D_B = bt_d_b->storage(); Real * D = tangent_stiffness_matrix->storage(); UInt offset_bt_d_b = bt_d_b_size * bt_d_b_size; UInt offset_d = tangent_size * tangent_size; for (UInt e = 0; e < nb_element; ++e) { Real * shapes_derivatives_val = shapes_derivatives.values + elem_filter_val[e]*size_of_shapes_derivatives*nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q) { transferBMatrixToSymVoigtBMatrix(shapes_derivatives_val, B, nb_nodes_per_element); Math::matrixt_matrix(bt_d_b_size, tangent_size, tangent_size, B, D, Bt_D); Math::matrix_matrix(bt_d_b_size, bt_d_b_size, tangent_size, Bt_D, B, Bt_D_B); shapes_derivatives_val += size_of_shapes_derivatives; D += offset_d; Bt_D_B += offset_bt_d_b; } } delete [] B; delete [] Bt_D; delete tangent_stiffness_matrix; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Vector * K_e = new Vector(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(); } /* -------------------------------------------------------------------------- */ 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; 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(); } /* -------------------------------------------------------------------------- */ 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::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); Vector nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += model->getDisplacement(); for(UInt gt = (UInt) _not_ghost; gt < (UInt) _casper; ++gt) { GhostType ghost_type = (GhostType) gt; 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 Vector & 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; Vector & 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::initInterpolateElementalField() { + AKANTU_DEBUG_IN(); + + const Mesh & mesh = model->getFEM().getMesh(); + + const Vector & position = mesh.getNodes(); + + Mesh::type_iterator it = mesh.firstType(spatial_dimension); + Mesh::type_iterator last = mesh.lastType(spatial_dimension); + + for (; it != last; ++it) { + UInt nb_element = mesh.getNbElement(*it); + if (nb_element == 0) continue; + + /// compute quadrature point position + UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(*it); + UInt nb_quad = nb_quad_per_element * nb_element; + Vector quad_coordinates(nb_quad, spatial_dimension); + + model->getFEM().interpolateOnQuadraturePoints(position, + quad_coordinates, + spatial_dimension, + *it); + + interpolationInvCoordinates.alloc(nb_quad, nb_quad_per_element, *it); + + Vector & interp_inv_coord = interpolationInvCoordinates(*it); + + ElementType type = *it; + +#define INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ + initInterpolation(quad_coordinates, interp_inv_coord) \ + + AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_INTERPOLATE_ELEMENTAL_FIELD); +#undef INIT_INTERPOLATE_ELEMENTAL_FIELD + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void Material::initInterpolation(const Vector & quad_coordinates, + Vector & interp_inv_coord) { + AKANTU_DEBUG_IN(); + + Vector & elem_fil = element_filter(type); + UInt nb_element = elem_fil.getSize(); + UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); + + types::Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); + + UInt quad_block = nb_quad_per_element * spatial_dimension; + UInt nb_quad2 = nb_quad_per_element * nb_quad_per_element; + + /// loop over the elements of the current material and element type + for (UInt el = 0; el < nb_element; ++el) { + + /// create a matrix containing the quadrature points coordinates + types::Matrix quad_coords(quad_coordinates.storage() + + quad_block * elem_fil(el), + nb_quad_per_element, + spatial_dimension); + + /// insert the quad coordinates in a matrix compatible with the interpolation + buildInterpolationCoodinates(quad_coords, quad_coord_matrix); + + /// create a matrix to store the matrix inversion result + types::Matrix inv_quad_coord_matrix(interp_inv_coord.storage() + + nb_quad2 * el, + nb_quad_per_element, + nb_quad_per_element); + + /// invert the interpolation matrix + Math::inv(nb_quad_per_element, + quad_coord_matrix.storage(), + inv_quad_coord_matrix.storage()); + } + + + AKANTU_DEBUG_OUT(); +} + /* -------------------------------------------------------------------------- */ void Material::interpolateStress(const ElementType type, - const Vector & quad_coordinates, const Vector & coordinates, Vector & result) { AKANTU_DEBUG_IN(); #define INTERPOLATE_ELEMENTAL_FIELD(type) \ interpolateElementalField(stress(type), \ - quad_coordinates, \ coordinates, \ result) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD); -#undef INTERPOLATE_ELEMENTAL_FIELD +#undef INTERPOLATE_ELEMENTAL_FIELD AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::interpolateElementalField(const Vector & field, - const Vector & quad_coordinates, const Vector & coordinates, Vector & result) { AKANTU_DEBUG_IN(); - UInt nb_element = element_filter(type).getSize(); + Vector & elem_fil = element_filter(type); + UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); - types::Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); - types::Matrix inv_quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); types::Matrix coefficients(nb_quad_per_element, field.getNbComponent()); + const Vector & interp_inv_coord = interpolationInvCoordinates(type); + Vector::const_iterator field_it = field.begin_reinterpret(nb_quad_per_element, field.getNbComponent(), nb_element, nb_quad_per_element * field.getNbComponent()); AKANTU_DEBUG_ASSERT(coordinates.getSize() % nb_element == 0, "Can't interpolate elemental field on elements, the coordinates vector has a wrong size"); UInt nb_points_per_elem = coordinates.getSize() / nb_element; types::Matrix coord_matrix(nb_points_per_elem, nb_quad_per_element); Vector::iterator result_it = result.begin_reinterpret(nb_points_per_elem, field.getNbComponent(), nb_element, nb_points_per_elem * field.getNbComponent()); - UInt quad_block = nb_quad_per_element * spatial_dimension; UInt coord_block = nb_points_per_elem * spatial_dimension; + UInt nb_quad2 = nb_quad_per_element * nb_quad_per_element; /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++result_it) { - /// create a matrix containing the quadrature points coordinates - types::Matrix quad_coords(quad_coordinates.storage() - + quad_block * element_filter(type)(el), - nb_quad_per_element, - spatial_dimension); - - /// insert the quad coordinates in a matrix compatible with the interpolation - buildInterpolationCoodinates(quad_coords, quad_coord_matrix); + /** + * create a matrix containing the inversion of the quadrature + * points' coordinates + * + */ + types::Matrix inv_quad_coord_matrix(interp_inv_coord.storage() + + nb_quad2 * el, + nb_quad_per_element, + nb_quad_per_element); - /// invert the interpolation matrix - Math::inv(nb_quad_per_element, - quad_coord_matrix.storage(), - inv_quad_coord_matrix.storage()); - /// multiply it by the field values over quadrature points to get - /// the interpolation coefficients + /** + * multiply it by the field values over quadrature points to get + * the interpolation coefficients + * + */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// create a matrix containing the points' coordinates types::Matrix coord(coordinates.storage() - + coord_block * element_filter(type)(el), - nb_points_per_elem, - spatial_dimension); + + coord_block * elem_fil(el), + nb_points_per_elem, + spatial_dimension); /// insert the points coordinates in a matrix compatible with the interpolation buildInterpolationCoodinates(coord, coord_matrix); /// multiply the coordinates matrix by the coefficients matrix and store the result (*result_it).mul(coord_matrix, coefficients); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const Vector & Material::getVector(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 << id << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getVector(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" <(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index b2820b750..6119a2537 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,460 +1,470 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 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 "fem.hh" //#include "mesh.hh" #include "data_accessor.hh" //#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; class CommunicationBuffer; } __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, * Vector & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : protected Memory, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(Vector & displacement, GhostType ghost_type = _not_ghost); void assembleResidual(GhostType ghost_type); /// compute the residual for this material virtual void computeStress(Vector & current_position, GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix void assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(Real h, const Element & element = ElementNull) = 0; /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// 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 virtual void interpolateStress(const ElementType type, - const Vector & quad_coordinates, const Vector & coordinates, Vector & result); + /// function to initialize the elemental field interpolation + /// function by inverting the quadrature points' coordinates + virtual void initInterpolateElementalField(); + protected: /// constitutive law virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(ElementType el_type, GhostType ghost_type = _not_ghost) {}; // /// constitutive law // virtual void computeNonLocalStress(ElementType el_type, // GhostType ghost_type = _not_ghost) { // AKANTU_DEBUG_TO_IMPLEMENT(); // }; /// compute the tangent stiffness matrix virtual void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type = _not_ghost) = 0; /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); template void assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const; /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Vector & field, - const Vector & quad_coordinates, const Vector & coordinates, Vector & result); + /// template function to initialize the elemental field interpolation + template + void initInterpolation(const Vector & quad_coordinates, + Vector & interp_inv_coord); + /// build the coordinate matrix for the interpolation on elemental field template inline void buildInterpolationCoodinates(const types::Matrix & coordinates, types::Matrix & coordMatrix); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for on element inline void computePotentialEnergyOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & epot); public: /// allocate an internal vector template void initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind = _ek_regular) const; /// resize an internal vector template void resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind = _ek_regular) const; /* ------------------------------------------------------------------------ */ template inline void gradUToF(const types::Matrix & grad_u, types::Matrix & F); inline void rightCauchy(const types::Matrix & F, types::Matrix & C); inline void leftCauchy (const types::Matrix & F, types::Matrix & B); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataToPack(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline UInt getNbDataToUnpack(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToPack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToUnpack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { } virtual void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) const { } virtual inline void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) { } virtual void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) { } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); virtual Real getEnergy(std::string type); 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); const Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; virtual Real getParam(const ID & param) const; virtual void setParam(const ID & param, Real value); protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material ID id; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// is the vector for potential energy initialized // bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// tell if using in non local mode or not bool is_non_local; /// spatial dimension UInt spatial_dimension; + /// elemental field interpolation coefficients + ByElementTypeReal interpolationInvCoordinates; + private: /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "material_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Vector::iterator stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Vector::iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++tangent_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ #define AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST \ ((stress_wf, StressBasedWeightFunction)) \ ((damage_wf, DamageWeightFunction )) \ ((base_wf, BaseWeightFunction )) #define AKANTU_MATERIAL_LIST \ ((2, (elastic , MaterialElastic ))) \ ((2, (viscoelastic , MaterialViscoElastic ))) \ ((2, (elastic_orthotropic , MaterialElasticOrthotropic ))) \ ((2, (elastic_caughey , MaterialElasticCaughey ))) \ ((2, (neohookean , MaterialNeohookean ))) \ ((2, (damage_linear , MaterialDamageLinear ))) \ ((2, (marigo , MaterialMarigo ))) \ ((2, (mazars , MaterialMazars ))) \ ((2, (vreepeerlings , MaterialVreePeerlings ))) \ ((2, (marigo_non_local , MaterialMarigoNonLocal ))) \ ((2, (mazars_non_local , MaterialMazarsNonLocal ))) \ ((2, (vreepeerlings_non_local, MaterialVreePeerlingsNonLocal))) \ ((2, (cohesive_bilinear , MaterialCohesiveBilinear ))) \ ((2, (cohesive_linear , MaterialCohesiveLinear ))) \ ((2, (cohesive_linear_extrinsic, MaterialCohesiveLinearExtrinsic ))) \ ((2, (cohesive_linear_exponential_extrinsic, MaterialCohesiveLinearExponentialExtrinsic ))) // ((3, (marigo_non_local_giry , MaterialMarigoNonLocal, (StressBasedWeigthFunction)))) #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #if defined(__INTEL_COMPILER) #pragma warning ( push ) /* warning #654: overloaded virtual function "akantu::Material::computeStress" is only partially overridden in class "akantu::Material*" */ #pragma warning ( disable : 654 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ // elastic materials #include "material_elastic.hh" #include "material_elastic_caughey.hh" #include "material_viscoelastic.hh" #include "material_neohookean.hh" #include "material_elastic_orthotropic.hh" // damage materials #include "material_damage.hh" #include "material_marigo.hh" #include "material_mazars.hh" #include "material_damage_linear.hh" #include "material_vreepeerlings.hh" #include "material_marigo_non_local.hh" #include "material_mazars_non_local.hh" #include "material_vreepeerlings_non_local.hh" // cohesive materials #include "material_cohesive.hh" #include "material_cohesive_linear.hh" #include "material_cohesive_bilinear.hh" #include "material_cohesive_linear_extrinsic.hh" #include "material_cohesive_linear_exponential_extrinsic.hh" #if defined(__INTEL_COMPILER) #pragma warning ( pop ) #endif //defined(__INTEL_COMPILER) #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc index fd7308c1e..6c8b40b46 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc @@ -1,280 +1,280 @@ /** * @file material_cohesive_linear.cc * @author Marco Vocialta * @date Mon Feb 20 12:14:16 2012 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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_cohesive_linear_extrinsic.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearExtrinsic::MaterialCohesiveLinearExtrinsic(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id), sigma_c_eff("sigma_c_eff",id), - delta_c("delta_c",id) { + delta_c("delta_c",id), + normal_stress(spatial_dimension), + tangential_stress(spatial_dimension) { AKANTU_DEBUG_IN(); sigma_c = 0; beta = 0; G_cI = 0; G_cII = 0; rand = 0; initInternalVector(sigma_c_eff, 1, _ek_cohesive); initInternalVector(delta_c, 1, _ek_cohesive); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearExtrinsic::~MaterialCohesiveLinearExtrinsic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); kappa = G_cII / G_cI; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::resizeCohesiveVectors() { MaterialCohesive::resizeCohesiveVectors(); resizeInternalVector(sigma_c_eff, _ek_cohesive); resizeInternalVector(delta_c, _ek_cohesive); FEM & fem_cohesive = model->getFEM("CohesiveFEM"); const Mesh & mesh = fem_cohesive.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { const Vector & elem_filter = element_filter(*it, _not_ghost); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = fem_cohesive.getNbQuadraturePoints(*it, _not_ghost); UInt nb_element_old = nb_element - sigma_insertion.getSize() / nb_quadrature_points; Vector & sigma_c_eff_vec = sigma_c_eff(*it, _not_ghost); Vector & delta_c_vec = delta_c(*it, _not_ghost); for (UInt el = nb_element_old; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Real new_sigma = sigma_insertion((el - nb_element_old) * nb_quadrature_points+q); Real new_delta = 2 * G_cI / new_sigma; sigma_c_eff_vec(el * nb_quadrature_points + q) = new_sigma; delta_c_vec(el * nb_quadrature_points + q) = new_delta; } } } } /* -------------------------------------------------------------------------- */ template bool MaterialCohesiveLinearExtrinsic::setParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); if(key == "sigma_c") { sstr >> sigma_c; } else if(key == "beta") { sstr >> beta; } else if(key == "G_cI") { sstr >> G_cI; } else if(key == "G_cII") { sstr >> G_cII; } else if(key == "rand") { sstr >> rand; } else { return Material::setParam(key, value, id); } return true; } /* -------------------------------------------------------------------------- */ template Real MaterialCohesiveLinearExtrinsic::computeEffectiveNorm(const types::Matrix & stress, const types::RVector & normal, const types::RVector & tangent) { AKANTU_DEBUG_IN(); Real normal_contrib, tangent_contrib; - types::RVector normal_stress(spatial_dimension); - types::RVector tangential_stress(spatial_dimension); normal_stress.mul(stress, normal); tangential_stress.mul(stress, tangent); normal_contrib = normal_stress.dot(normal); tangent_contrib = tangential_stress.dot(tangent); if (normal_contrib < 0) normal_contrib = 0; AKANTU_DEBUG_OUT(); return std::sqrt(normal_contrib*normal_contrib + tangent_contrib*tangent_contrib/beta/beta); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Vector::iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Vector::iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Vector::const_iterator normal_it = normal.begin(spatial_dimension); Vector::iteratortraction_end = tractions(el_type, ghost_type).end(spatial_dimension); Vector::iteratorsigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Vector::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); Vector::iteratordelta_c_it = delta_c(el_type, ghost_type).begin(); Vector::iteratordamage_it = damage(el_type, ghost_type).begin(); /// compute scalars Real beta2_kappa2 = beta*beta/kappa/kappa; Real beta2_kappa = beta*beta/kappa; /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); types::Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; types::Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2_kappa2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); /// full damage case or zero displacement case if (delta >= *delta_c_it || delta == 0) { /// set traction to zero (*traction_it).clear(); *damage_it = delta >= *delta_c_it; *delta_max_it = *damage_it * (*delta_c_it); } /// element not fully damaged else { /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ *traction_it = tangential_opening; *traction_it *= beta2_kappa; *traction_it += normal_opening; /// update maximum displacement *delta_max_it = std::max(*delta_max_it, delta); *damage_it = *delta_max_it / *delta_c_it; Real k = *sigma_c_it / *delta_max_it * (1. - *damage_it); *traction_it *= k; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_cohesive_linear> [" << std::endl; stream << space << " + sigma_c : " << sigma_c << std::endl; stream << space << " + beta : " << beta << std::endl; stream << space << " + G_cI : " << G_cI << std::endl; stream << space << " + G_cII : " << G_cII << std::endl; stream << space << " + rand : " << rand << std::endl; if(this->isInit()) { stream << space << " + kappa : " << kappa << std::endl; } MaterialCohesive::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveLinearExtrinsic); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh index da3c298ab..4167beb96 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh @@ -1,131 +1,137 @@ /** * @file material_cohesive_linear.hh * @author Marco Vocialta * @date Mon Feb 20 12:00:34 2012 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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_cohesive.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material linear damage for extrinsic case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: 0) * - G_cI : fracture energy for mode I (default: 0) * - G_cII : fracture energy for mode II (default: 0) * - rand : randomness factor (default: 0) */ template class MaterialCohesiveLinearExtrinsic : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveLinearExtrinsic(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialCohesiveLinearExtrinsic(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set patameters virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// initialize the material computed parameter virtual void initMaterial(); /// resize vectors for new cohesive elements virtual void resizeCohesiveVectors(); /// compute effective stress norm for insertion check virtual Real computeEffectiveNorm(const types::Matrix & stress, const types::RVector & normal, const types::RVector & tangent); protected: /// constitutive law void computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// critical effective stress ByElementTypeReal sigma_c_eff; /// beta parameter Real beta; /// mode I fracture energy Real G_cI; /// mode II fracture energy Real G_cII; /// kappa parameter Real kappa; /// critical displacement ByElementTypeReal delta_c; + /// vector to temporarily store the normal stress for the norm + types::RVector normal_stress; + + /// vector to temporarily store the tangential stress for the norm + types::RVector tangential_stress; + }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_cohesive_linear_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc index 12b89acc8..6a5df6ad3 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc @@ -1,932 +1,922 @@ /** * @file solid_mechanics_model_cohesive.cc * @author Marco Vocialta * @date Thu Apr 19 10:42:11 2012 * * @brief Solid mechanics model for cohesive elements * * @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 #include #include #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), mesh_facets(mesh.getSpatialDimension(), mesh.getNodes().getID(), id, memory_id), - quad_elements("quad_elements", id), elements_quad_facets("elements_quad_facets", id), facet_stress(0, spatial_dimension * spatial_dimension, "facet_stress"), facets_to_cohesive_el(0, 2, "facets_to_cohesive_el"), fragment_to_element("fragment_to_element", id) { AKANTU_DEBUG_IN(); registerFEMObject("CohesiveFEM", mesh, spatial_dimension); MeshUtils::buildAllFacets(mesh, mesh_facets); /// assign cohesive type Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { const Vector & connectivity = mesh_facets.getConnectivity(*it); if (connectivity.getSize() != 0) { type_facet = *it; break; } } if (type_facet == _segment_2) type_cohesive = _cohesive_2d_4; else if (type_facet == _segment_3) type_cohesive = _cohesive_2d_6; mesh.addConnectivityType(type_cohesive); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFull(std::string material_file, AnalysisMethod method) { SolidMechanicsModel::initFull(material_file, method); if(material_file != "") { initCohesiveMaterial(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initCohesiveMaterial() { AKANTU_DEBUG_IN(); /// find the cohesive index cohesive_index = 0; while ((dynamic_cast(materials[cohesive_index]) == NULL) && cohesive_index <= materials.size()) ++cohesive_index; AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(), "No cohesive materials in the material input file"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModelCohesive::initModel() { SolidMechanicsModel::initModel(); getFEM("CohesiveFEM").initShapeFunctions(_not_ghost); getFEM("CohesiveFEM").initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initExtrinsic() { AKANTU_DEBUG_IN(); /// assign sigma_c to each facet MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); UInt nb_facet = mesh_facets.getNbElement(type_facet); sigma_lim.resize(nb_facet); const Real sigma_c = mat_cohesive->getSigmaC(); const Real rand = mat_cohesive->getRandFactor(); std::srand(time(NULL)); if (facets_check.getSize() < 1) { const Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); facets_check.resize(nb_facet); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)(1) != ElementNull) { facets_check(f) = true; sigma_lim(f) = sigma_c * (1 + std::rand()/(Real)RAND_MAX * rand); } else { facets_check(f) = false; sigma_lim(f) = std::numeric_limits::max(); } } } /// compute normals on facets registerFEMObject("FacetsFEM", mesh_facets, spatial_dimension-1); getFEM("FacetsFEM").initShapeFunctions(); getFEM("FacetsFEM").computeNormalsOnControlPoints(); /// THIS HAS TO BE CHANGED: /* ------------------------------------------------------------------------ */ const Vector & normals = getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet); Vector::const_iterator normal_it = normals.begin(spatial_dimension); tangents.resize(normals.getSize()); tangents.extendComponentsInterlaced(spatial_dimension, tangents.getNbComponent()); Vector::iterator tangent_it = tangents.begin(spatial_dimension); for (UInt i = 0; i < normals.getSize(); ++i, ++normal_it, ++tangent_it) { Math::normal2( (*normal_it).storage(), (*tangent_it).storage() ); } /* ------------------------------------------------------------------------ */ /// compute quadrature points coordinates on facets Vector & position = mesh.getNodes(); UInt nb_quad_per_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_tot_quad = nb_quad_per_facet * nb_facet; Vector quad_facets(nb_tot_quad, spatial_dimension); getFEM("FacetsFEM").interpolateOnQuadraturePoints(position, quad_facets, spatial_dimension, type_facet); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; - /// compute quadrature point position - UInt nb_quad_per_element = getFEM().getNbQuadraturePoints(*it); - UInt nb_quad = nb_quad_per_element * nb_element; - quad_elements.alloc(nb_quad, spatial_dimension, *it); - - getFEM().interpolateOnQuadraturePoints(position, - quad_elements(*it), - spatial_dimension, - *it); - /// compute elements' quadrature points and list of facet /// quadrature points positions by element Vector & facet_to_element = mesh_facets.getSubelementToElement(*it); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); elements_quad_facets.alloc(nb_element * nb_facet_per_elem * nb_quad_per_facet, spatial_dimension, *it); Vector & el_q_facet = elements_quad_facets(*it); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { UInt global_facet = facet_to_element(el, f).element; for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_facets(global_facet * nb_quad_per_facet + q, s); } } } } } + /// initialize the interpolation function + materials[0]->initInterpolateElementalField(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); Vector facet_insertion; UInt nb_materials = materials.size(); UInt nb_quad_per_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = mesh_facets.getNbElement(type_facet); /// list of stresses on facet quadrature points for every element ByElementTypeReal stress_on_facet; /// vector containing stresses coming from the two elements of each facets facet_stress.resize(2 * nb_facet * nb_quad_per_facet); Vector facet_stress_count(nb_facet); facet_stress_count.clear(); /// loop over materials for (UInt m = 0; m < nb_materials; ++m) { if (m == cohesive_index) continue; Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); /// loop over element type for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; - const Vector & quad_elem = quad_elements(*it); const Vector & el_q_facet = elements_quad_facets(*it); stress_on_facet.alloc(el_q_facet.getSize(), spatial_dimension * spatial_dimension, *it); Vector & stress_on_f = stress_on_facet(*it); /// interpolate stress on facet quadrature points positions materials[m]->interpolateStress(*it, - quad_elem, el_q_facet, stress_on_f); /// store the interpolated stresses on the facet_stress vector Vector & facet_to_element = mesh_facets.getSubelementToElement(*it); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { UInt global_facet = facet_to_element(el, f).element; for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { facet_stress(2 * global_facet * nb_quad_per_facet + facet_stress_count(global_facet) * nb_quad_per_facet + q, s) = stress_on_f(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s); } } facet_stress_count(global_facet) = true; } } } } MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); mat_cohesive->checkInsertion(facet_stress, facet_insertion); if (facet_insertion.getSize() != 0) insertCohesiveElements(facet_insertion); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertCohesiveElements(const Vector & facet_insertion) { AKANTU_DEBUG_IN(); if(facet_insertion.getSize() == 0) return; Vector doubled_nodes(0, 2); Vector doubled_facets(0, 2); /// update mesh for (UInt f = 0; f < facet_insertion.getSize(); ++f) { Element facet(type_facet, facet_insertion(f)); doubleFacet(facet, doubled_nodes, doubled_facets); } /// double middle nodes if it's the case if (type_facet == _segment_3) doubleMiddleNode(doubled_nodes, doubled_facets); /// loop over doubled facets to insert cohesive elements Vector & conn_cohesive = mesh.getConnectivity(type_cohesive); const Vector & conn_facet = mesh_facets.getConnectivity(type_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); const Vector & position = mesh.getNodes(); Vector & element_mat = getElementMaterial(type_cohesive); Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); for (UInt f = 0; f < doubled_facets.getSize(); ++f) { UInt nb_cohesive_elements = conn_cohesive.getSize(); conn_cohesive.resize(nb_cohesive_elements + 1); UInt first_facet = doubled_facets(f, 0); UInt second_facet = doubled_facets(f, 1); /// copy first facet's connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) conn_cohesive(nb_cohesive_elements, n) = conn_facet(first_facet, n); /// check if first nodes of the two facets are coincident or not UInt first_facet_node = conn_facet(first_facet, 0); UInt second_facet_node = conn_facet(second_facet, 0); bool second_facet_inversion = false; for (UInt dim = 0; dim < mesh.getSpatialDimension(); ++dim) { if (position(first_facet_node, dim) != position(second_facet_node, dim)) { second_facet_inversion = true; break; } } /// if the two nodes are coincident second facet connectivity is /// normally copied, otherwise the copy is reverted if (!second_facet_inversion) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet) = conn_facet(second_facet, n); } else { for (UInt n = 0; n < nb_nodes_per_facet; ++n) conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet) = conn_facet(second_facet, nb_nodes_per_facet - n - 1); } /// assign the cohesive material to the new element element_mat.resize(nb_cohesive_elements + 1); element_mat(nb_cohesive_elements) = cohesive_index; materials[cohesive_index]->addElement(type_cohesive, nb_cohesive_elements, _not_ghost); /// update element_to_facet vectors Element cohesive_element(type_cohesive, nb_cohesive_elements, _not_ghost, _ek_cohesive); element_to_facet(first_facet)(1) = cohesive_element; element_to_facet(second_facet)(1) = cohesive_element; /// update facets_to_cohesive_el vector facets_to_cohesive_el.resize(nb_cohesive_elements + 1); facets_to_cohesive_el(nb_cohesive_elements, 0) = first_facet; facets_to_cohesive_el(nb_cohesive_elements, 1) = second_facet; } /// update shape functions getFEM("CohesiveFEM").initShapeFunctions(_not_ghost); /// resize cohesive material vectors MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); AKANTU_DEBUG_ASSERT(mat_cohesive, "No cohesive materials in the materials vector"); mat_cohesive->resizeCohesiveVectors(); /// update nodal values updateDoubledNodes(doubled_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::doubleMiddleNode(Vector & doubled_nodes, const Vector & doubled_facets) { AKANTU_DEBUG_IN(); Vector & conn_facet = mesh_facets.getConnectivity(type_facet); Vector & position = mesh.getNodes(); Vector > & elem_to_facet = mesh_facets.getElementToSubelement(type_facet); for (UInt f = 0; f < doubled_facets.getSize(); ++f) { UInt facet_first = doubled_facets(f, 0); UInt facet_second = doubled_facets(f, 1); UInt new_node = position.getSize(); /// store doubled nodes UInt nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(nb_doubled_nodes + 1); UInt old_node = conn_facet(facet_first, 2); doubled_nodes(nb_doubled_nodes, 0) = old_node; doubled_nodes(nb_doubled_nodes, 1) = new_node; /// update position position.resize(new_node + 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) position(new_node, dim) = position(old_node, dim); /// update facet connectivity conn_facet(facet_second, 2) = new_node; /// update element connectivity for (UInt el = 0; el < elem_to_facet(facet_second).getSize(); ++el) { const ElementType type_elem = elem_to_facet(facet_second)(el).type; if (type_elem != _not_defined) { UInt elem_global = elem_to_facet(facet_second)(el).element; Vector & conn_elem = mesh.getConnectivity(type_elem); for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) { if (conn_elem(elem_global, n) == old_node) conn_elem(elem_global, n) = new_node; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateDoubledNodes(const Vector & doubled_nodes) { AKANTU_DEBUG_IN(); UInt nb_old_nodes = displacement->getSize(); UInt nb_new_nodes = nb_old_nodes + doubled_nodes.getSize(); displacement ->resize(nb_new_nodes); velocity ->resize(nb_new_nodes); acceleration ->resize(nb_new_nodes); increment_acceleration->resize(nb_new_nodes); force ->resize(nb_new_nodes); residual ->resize(nb_new_nodes); boundary ->resize(nb_new_nodes); mass ->resize(nb_new_nodes); /** * @todo temporary patch, should be done in a better way that works * always (pbc, parallel, ...) **/ delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); for (UInt n = 0; n < doubled_nodes.getSize(); ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); for (UInt dim = 0; dim < displacement->getNbComponent(); ++dim) { (*displacement)(new_node, dim) = (*displacement)(old_node, dim); } for (UInt dim = 0; dim < velocity->getNbComponent(); ++dim) { (*velocity)(new_node, dim) = (*velocity)(old_node, dim); } for (UInt dim = 0; dim < acceleration->getNbComponent(); ++dim) { (*acceleration)(new_node, dim) = (*acceleration)(old_node, dim); } for (UInt dim = 0; dim < increment_acceleration->getNbComponent(); ++dim) { (*increment_acceleration)(new_node, dim) = (*increment_acceleration)(old_node, dim); } } assembleMassLumped(); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::doubleFacet(Element & facet, Vector & doubled_nodes, Vector & doubled_facets) { AKANTU_DEBUG_IN(); const UInt f_index = facet.element; const ElementType type_facet = facet.type; const ElementType type_subfacet = mesh.getFacetElementType(type_facet); const UInt nb_subfacet = mesh.getNbFacetsPerElement(type_facet); Vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet); Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); Vector & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet); /// adding a new facet by copying original one /// create new connectivity Vector & conn_facet = mesh_facets.getConnectivity(type_facet); UInt nb_facet = conn_facet.getSize(); conn_facet.resize(nb_facet + 1); for (UInt n = 0; n < conn_facet.getNbComponent(); ++n) conn_facet(nb_facet, n) = conn_facet(f_index, n); /// store doubled facets UInt nb_doubled_facets = doubled_facets.getSize(); doubled_facets.resize(nb_doubled_facets + 1); doubled_facets(nb_doubled_facets, 0) = f_index; doubled_facets(nb_doubled_facets, 1) = nb_facet; /// update elements connected to facet Vector first_facet_list = element_to_facet(f_index); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(f_index)(1) = ElementNull; element_to_facet(nb_facet)(0) = element_to_facet(nb_facet)(1); element_to_facet(nb_facet)(1) = ElementNull; /// update facet_to_element vector ElementType type = element_to_facet(nb_facet)(0).type; UInt el = element_to_facet(nb_facet)(0).element; Vector & facet_to_element = mesh_facets.getSubelementToElement(type); UInt i; for (i = 0; facet_to_element(el, i).element != f_index && i <= facet_to_element.getNbComponent(); ++i); facet_to_element(el, i).element = nb_facet; /// create new links to subfacets and update list of facets /// connected to subfacets subfacet_to_facet.resize(nb_facet + 1); for (UInt sf = 0; sf < nb_subfacet; ++sf) { subfacet_to_facet(nb_facet, sf) = subfacet_to_facet(f_index, sf); UInt sf_index = subfacet_to_facet(f_index, sf).element; /// find index to start looping around facets connected to current /// subfacet UInt start = 0; UInt nb_connected_facets = facet_to_subfacet(sf_index).getSize(); while (facet_to_subfacet(sf_index)(start).element != f_index && start <= facet_to_subfacet(sf_index).getSize()) ++start; /// add the new facet to the list next to the original one ++nb_connected_facets; facet_to_subfacet(sf_index).resize(nb_connected_facets); for (UInt f = nb_connected_facets - 1; f > start; --f) { facet_to_subfacet(sf_index)(f) = facet_to_subfacet(sf_index)(f - 1); } /// check if the new facet should be inserted before or after the /// original one: the second element connected to both original /// and new facet will be _not_defined, so I check if the first /// one is equal to one of the elements connected to the following /// facet in the facet_to_subfacet vector UInt f_start = facet_to_subfacet(sf_index)(start).element; UInt f_next; if (start + 2 == nb_connected_facets) f_next = facet_to_subfacet(sf_index)(0).element; else f_next = facet_to_subfacet(sf_index)(start + 2).element; if ((element_to_facet(f_start)(0) == element_to_facet(f_next)(0)) || ( element_to_facet(f_start)(0) == element_to_facet(f_next)(1))) facet_to_subfacet(sf_index)(start).element = nb_facet; else facet_to_subfacet(sf_index)(start + 1).element = nb_facet; /// loop on every facet connected to the current subfacet for (UInt f = start + 2; ; ++f) { /// reset f in order to continue looping from the beginning if (f == nb_connected_facets) f = 0; /// exit loop if it reaches the end if (f == start) break; /// if current loop facet is on the boundary, double subfacet UInt f_global = facet_to_subfacet(sf_index)(f).element; if (element_to_facet(f_global)(1).type == _not_defined || element_to_facet(f_global)(1).kind == _ek_cohesive) { doubleSubfacet(subfacet_to_facet(f_index, sf), start, f, doubled_nodes); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::doubleSubfacet(const Element & subfacet, UInt start, UInt end, Vector & doubled_nodes) { AKANTU_DEBUG_IN(); const UInt sf_index = subfacet.element; const ElementType type_subfacet = subfacet.type; Vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet); UInt nb_subfacet = facet_to_subfacet.getSize(); Vector & conn_point = mesh_facets.getConnectivity(_point); Vector & position = mesh.getNodes(); /// add the new subfacet if (spatial_dimension == 2) { UInt new_node = position.getSize(); /// add new node in connectivity UInt new_subfacet = conn_point.getSize(); conn_point.resize(new_subfacet + 1); conn_point(new_subfacet) = new_node; /// store doubled nodes UInt nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(nb_doubled_nodes + 1); UInt old_node = doubled_nodes(nb_doubled_nodes, 0) = conn_point(sf_index); doubled_nodes(nb_doubled_nodes, 1) = new_node; /// update position position.resize(new_node + 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) position(new_node, dim) = position(old_node, dim); } /// create a vector for the new subfacet in facet_to_subfacet facet_to_subfacet.resize(nb_subfacet + 1); UInt nb_connected_facets = facet_to_subfacet(sf_index).getSize(); /// loop over facets from start to end for (UInt f = start + 1; ; ++f) { /// reset f in order to continue looping from the beginning if (f == nb_connected_facets) f = 0; UInt f_global = facet_to_subfacet(sf_index)(f).element; ElementType type_facet = facet_to_subfacet(sf_index)(f).type; Vector & conn_facet = mesh_facets.getConnectivity(type_facet); UInt old_node = conn_point(sf_index); UInt new_node = conn_point(conn_point.getSize() - 1); /// update facet connectivity UInt i; for (i = 0; conn_facet(f_global, i) != old_node && i <= conn_facet.getNbComponent(); ++i); conn_facet(f_global, i) = new_node; /// update element connectivity Vector > & elem_to_facet = mesh_facets.getElementToSubelement(type_facet); for (UInt el = 0; el < elem_to_facet(f_global).getSize(); ++el) { const ElementType type_elem = elem_to_facet(f_global)(el).type; if (type_elem != _not_defined) { UInt elem_global = elem_to_facet(f_global)(el).element; Vector & conn_elem = mesh.getConnectivity(type_elem); for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) { if (conn_elem(elem_global, n) == old_node) conn_elem(elem_global, n) = new_node; } } } /// update subfacet_to_facet vector Vector & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet); for (i = 0; subfacet_to_facet(f_global, i).element != sf_index && i <= subfacet_to_facet.getNbComponent(); ++i); subfacet_to_facet(f_global, i).element = nb_subfacet; /// add current facet to facet_to_subfacet last position Element current_facet(type_facet, f_global); facet_to_subfacet(nb_subfacet).push_back(current_facet); /// exit loop if it reaches the end if (f == end) break; } /// rearrange the facets connected to the original subfacet and /// compute the new number of facets connected to it if (end < start) { for (UInt f = 0; f < start - end; ++f) facet_to_subfacet(sf_index)(f) = facet_to_subfacet(sf_index)(f + end + 1); nb_connected_facets = start - end; } else { for (UInt f = 1; f < nb_connected_facets - end; ++f) facet_to_subfacet(sf_index)(start + f) = facet_to_subfacet(sf_index)(end + f); nb_connected_facets -= end - start; } /// resize list of facets of the original subfacet facet_to_subfacet(sf_index).resize(nb_connected_facets); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::buildFragmentsList() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); const UInt max = std::numeric_limits::max(); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element != 0) { /// initialize the list of checked elements and the list of /// elements to be checked fragment_to_element.alloc(nb_element, 1, *it); for (UInt el = 0; el < nb_element; ++el) fragment_to_element(*it)(el) = max; } } Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); nb_fragment = 0; it = mesh.firstType(spatial_dimension); MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); const Vector & damage = mat_cohesive->getDamage(type_cohesive); UInt nb_quad_cohesive = getFEM("CohesiveFEM").getNbQuadraturePoints(type_cohesive); for (; it != last; ++it) { Vector & checked_el = fragment_to_element(*it); UInt nb_element = checked_el.getSize(); if (nb_element != 0) { Vector & facet_to_element = mesh_facets.getSubelementToElement(*it); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); /// loop on elements for (UInt el = 0; el < nb_element; ++el) { if (checked_el(el) == max) { /// build fragment ++nb_fragment; checked_el(el) = nb_fragment - 1; Vector elem_to_check; Element current_el(*it, el); elem_to_check.push_back(current_el); /// keep looping while there are elements to check while (elem_to_check.getSize() != 0) { UInt nb_elem_check = elem_to_check.getSize(); for (UInt el_check = 0; el_check < nb_elem_check; ++el_check) { Element current_el = elem_to_check(el_check); for (UInt f = 0; f < nb_facet_per_elem; ++f) { /// find adjacent element on current facet UInt global_facet = facet_to_element(current_el.element, f).element; Element next_el; for (UInt i = 0; i < 2; ++i) { next_el = element_to_facet(global_facet)(i); if (next_el != current_el) break; } if (next_el.kind == _ek_cohesive) { /// fragmention occurs when the cohesive element has /// reached damage = 1 on every quadrature point UInt q = 0; while (q < nb_quad_cohesive && damage(next_el.element * nb_quad_cohesive + q) == 1) ++q; if (q == nb_quad_cohesive) next_el = ElementNull; else { /// check which facet is the correct one UInt other_facet_index = facets_to_cohesive_el(next_el.element, 0) == global_facet; UInt other_facet = facets_to_cohesive_el(next_el.element, other_facet_index); /// get the other regualar element next_el = element_to_facet(other_facet)(0); } } /// if it exists, add it to the fragment list if (next_el != ElementNull) { Vector & checked_next_el = fragment_to_element(next_el.type); /// check if the element isn't already part of a fragment if (checked_next_el(next_el.element) == max) { checked_next_el(next_el.element) = nb_fragment - 1; elem_to_check.push_back(next_el); } } } } /// erase elements that have already been checked for (UInt el_check = nb_elem_check; el_check > 0; --el_check) elem_to_check.erase(el_check - 1); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); return mat_cohesive->getReversibleEnergy(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModelCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); return mat_cohesive->getDissipatedEnergy(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh index c46503b86..b3ed0942e 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh @@ -1,210 +1,207 @@ /** * @file solid_mechanics_model_cohesive.hh * @author Marco Vocialta * @date Thu Apr 19 10:34:00 2012 * * @brief Solid mechanics model for cohesive elements * * @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 "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ __BEGIN_AKANTU__ class SolidMechanicsModelCohesive : public SolidMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate< IntegratorCohesive, ShapeCohesive > MyFEMCohesiveType; SolidMechanicsModelCohesive(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "solid_mechanics_model_cohesive", const MemoryID & memory_id = 0); // virtual ~SolidMechanicsModelCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function to perform a stress check on each facet and insert /// cohesive elements if needed void checkCohesiveStress(); /// function to insert cohesive elements on the selected facets void insertCohesiveElements(const Vector & facet_insertion); /// initialize completely the model void initFull(std::string material_file, AnalysisMethod method = _explicit_dynamic); /// initialize completely the model for extrinsic elements void initExtrinsic(); /// initialize the model void initModel(); /// initialize cohesive material void initCohesiveMaterial(); /// build fragments list void buildFragmentsList(); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); private: /// function to double a given facet and update the list of doubled /// nodes void doubleFacet(Element & facet, Vector & doubled_nodes, Vector & doubled_facets); /// function to double a subfacet given start and end index for /// local facet_to_subfacet vector, and update the list of doubled /// nodes void doubleSubfacet(const Element & subfacet, UInt start, UInt end, Vector & doubled_nodes); /// double middle nodes if facets are _segment_3 void doubleMiddleNode(Vector & doubled_nodes, const Vector & doubled_facets); /// function to update nodal parameters for doubled nodes void updateDoubledNodes(const Vector & doubled_nodes); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the cohesive element type AKANTU_GET_MACRO(CohesiveElementType, type_cohesive, ElementType); /// get the cohesive index AKANTU_GET_MACRO(CohesiveIndex, cohesive_index, UInt); /// get the facet type AKANTU_GET_MACRO(FacetType, type_facet, ElementType); /// get the facet mesh AKANTU_GET_MACRO(MeshFacets, mesh_facets, const Mesh &); /// get the sigma limit vector for automatic insertion AKANTU_GET_MACRO(SigmaLimit, sigma_lim, const Vector &); AKANTU_GET_MACRO_NOT_CONST(SigmaLimit, sigma_lim, Vector &); /// get the facets check vector AKANTU_GET_MACRO_NOT_CONST(FacetsCheck, facets_check, Vector &); /// get the stress on facets vector AKANTU_GET_MACRO(StressOnFacets, facet_stress, const Vector &); /// get the number of fragments AKANTU_GET_MACRO(NbFragment, nb_fragment, UInt); /// get the fragment_to_element vectors AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FragmentToElement, fragment_to_element, UInt); /// THIS HAS TO BE CHANGED AKANTU_GET_MACRO(Tangents, tangents, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// cohesive element type ElementType type_cohesive; /// facet type ElementType type_facet; /// cohesive material index in materials vector UInt cohesive_index; /// mesh containing facets and their data structures Mesh mesh_facets; /// vector containing a sigma limit for automatic insertion Vector sigma_lim; /// vector containing facets in which cohesive elements can be automatically inserted Vector facets_check; /// @todo store tangents when normals are computed: Vector tangents; - /// quadrature points coordinates by element - ByElementTypeReal quad_elements; - /// list of facet quadrature points positions by element ByElementTypeReal elements_quad_facets; /// stress on facets on the two sides by quadrature point Vector facet_stress; /// list of facets connected to each cohesive element Vector facets_to_cohesive_el; /// fragment number for each element ByElementTypeUInt fragment_to_element; /// number of fragments UInt nb_fragment; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "solid_mechanics_model_cohesive_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModelCohesive & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc index 8c1abe1fb..fbae84139 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc @@ -1,174 +1,175 @@ /** * @file test_interpolate_stress.cc * @author Marco Vocialta * @date Wed Jun 6 14:25:54 2012 * * @brief Test for the stress interpolation function * * @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 #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "io_helper.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; Real plane(Real x, Real y) { return 1. + 2. * x + 3. * y; } int main(int argc, char *argv[]) { initialize(argc, argv); debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 2; const ElementType type = _triangle_6; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("triangle.msh", mesh); const ElementType type_facet = mesh.getFacetElementType(type);; MeshUtils::buildAllFacets(mesh, mesh); SolidMechanicsModel model(mesh); /// model initialization model.initFull("../material.dat"); Vector & position = mesh.getNodes(); UInt nb_facet = mesh.getNbElement(type_facet); UInt nb_element = mesh.getNbElement(type); /// compute quadrature points positions on facets typedef FEMTemplate MyFEMType; model.registerFEMObject("FacetsFEM", mesh, spatial_dimension-1); model.getFEM("FacetsFEM").initShapeFunctions(); UInt nb_quad_per_facet = model.getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_tot_quad = nb_quad_per_facet * nb_facet; Vector quad_facets(nb_tot_quad, spatial_dimension); model.getFEM("FacetsFEM").interpolateOnQuadraturePoints(position, quad_facets, spatial_dimension, type_facet); Vector & facet_to_element = mesh.getSubelementToElement(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); Vector el_q_facet(nb_element * nb_facet_per_elem * nb_quad_per_facet, spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { UInt global_facet = facet_to_element(el, f).element; for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_facets(global_facet * nb_quad_per_facet + q, s); } } } } /// compute quadrature points position of the elements UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type); UInt nb_tot_quad_el = nb_quad_per_element * nb_element; Vector quad_elements(nb_tot_quad_el, spatial_dimension); model.getFEM().interpolateOnQuadraturePoints(position, quad_elements, spatial_dimension, type); /// assign some values to stresses Vector & stress = const_cast&>(model.getMaterial(0).getStress(type)); for (UInt q = 0; q < nb_tot_quad_el; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { stress(q, s) = s * plane(quad_elements(q, 0), quad_elements(q, 1)); } } /// interpolate stresses on facets' quadrature points + model.getMaterial(0).initInterpolateElementalField(); + Vector interpolated_stress(nb_element * nb_facet_per_elem * nb_quad_per_facet, stress.getNbComponent()); model.getMaterial(0).interpolateStress(type, - quad_elements, el_q_facet, interpolated_stress); /// check results for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { Real x = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 0); Real y = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 1); Real theoretical = s * plane(x, y); Real numerical = interpolated_stress(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s); Real tolerance = 1000 * (theoretical + 1) * std::numeric_limits::epsilon(); if (std::abs(theoretical - numerical) > tolerance) { std::cout << "Theoretical and numerical values aren't coincident!" << std::endl; return EXIT_FAILURE; } } } } } std::cout << "OK: Stress interpolation test passed." << std::endl; return EXIT_SUCCESS; }