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 new file mode 100644 index 000000000..b4439d991 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc @@ -0,0 +1,253 @@ +/** + * @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.hh" +#include "sparse_matrix.hh" +#include "dof_synchronizer.hh" + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +MaterialCohesiveLinearExtrinsic::MaterialCohesiveLinearExtrinsic(Model & model, const ID & id) : + MaterialCohesive(model,id), + sigma_c_eff("sigma_c_eff",id), + delta_max("delta max",id), + delta_c("delta_c",id) { + 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_max, 1, _ek_cohesive); + initInternalVector(delta_c, 1, _ek_cohesive); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +MaterialCohesiveLinearExtrinsic::~MaterialCohesiveLinearExtrinsic() { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialCohesiveLinearExtrinsic::initMaterial() { + AKANTU_DEBUG_IN(); + + MaterialCohesive::initMaterial(); + + kappa = G_cII / G_cI; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialCohesiveLinearExtrinsic::resizeCohesiveVectors(const Vector & sigma_insertion) { + MaterialCohesive::resizeCohesiveVectors(); + + resizeInternalVector(sigma_c_eff, _ek_cohesive); + resizeInternalVector(delta_max, _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_element_old = nb_element - sigma_insertion.getSize(); + UInt nb_quadrature_points = fem_cohesive.getNbQuadraturePoints(*it, _not_ghost); + + 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) { + Real new_sigma = sigma_insertion(el - nb_element_old); + Real new_delta = 2 * G_cI / new_sigma; + for (UInt q = 0; q < nb_quadrature_points; ++q) { + sigma_c_eff_vec(el * nb_quadrature_points + q) = new_sigma; + delta_c_vec(el * nb_quadrature_points + q) = new_delta; + } + } + } + +} + +/* -------------------------------------------------------------------------- */ +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; +} + +/* -------------------------------------------------------------------------- */ +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(); + + /// 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) { + + /// 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 + if (delta >= *delta_c_it || delta == 0) { + + /// set traction to zero + for (UInt i = 0; i < spatial_dimension; ++i) + (*traction_it)(i) = 0.; + + } + /// 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; + + /// crack opening case + if (delta > *delta_max_it) { + Real k = *sigma_c_it / delta * (1. - delta / *delta_c_it); + *traction_it *= k; + + /// update maximum displacement + *delta_max_it = delta; + + } + /// unloading-reloading case + else { + Real k = *sigma_c_it / *delta_max_it * (1. - *delta_max_it / *delta_c_it); + *traction_it *= k; + } + + } + + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +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(is_init) { + stream << space << " + kappa : " << kappa << std::endl; + } + MaterialCohesive::printself(stream, indent + 1); + stream << space << "]" << std::endl; +} +/* -------------------------------------------------------------------------- */ + + +__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 new file mode 100644 index 000000000..807a5b444 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh @@ -0,0 +1,149 @@ +/** + * @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) + */ + +class MaterialCohesiveLinearExtrinsic : public MaterialCohesive { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +public: + + MaterialCohesiveLinearExtrinsic(Model & 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(const Vector & sigma_insertion); + +protected: + + /// constitutive law + void computeTraction(const Vector & normal, + ElementType el_type, + GhostType ghost_type = _not_ghost); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + + /// get sigma_c + AKANTU_GET_MACRO(SigmaC, sigma_c, Real); + + /// get rand + AKANTU_GET_MACRO(RandFactor, rand, Real); + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +protected: + + /// critical stress + Real sigma_c; + + /// 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; + + /// maximum displacement + ByElementTypeReal delta_max; + + /// critical displacement + ByElementTypeReal delta_c; + + /// randomness factor + Real rand; + +}; + + +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ + +//#include "material_cohesive_linear_inline_impl.cc" + +/// standard output stream operator +inline std::ostream & operator <<(std::ostream & stream, const MaterialCohesiveLinearExtrinsic & _this) +{ + _this.printself(stream); + return stream; +} + + +__END_AKANTU__ + +#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc new file mode 100644 index 000000000..0a74bfb47 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc @@ -0,0 +1,194 @@ +/** + * @file material_elastic_orthotropic.cc + * @author Marco Vocialta + * @date Thu Apr 12 11:28:23 2012 + * + * @brief Orthotropic elastic material + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +#include "material_elastic_orthotropic.hh" +#include "solid_mechanics_model.hh" +#include + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +MaterialElasticOrthotropic::MaterialElasticOrthotropic(Model & model, const ID & id) : + Material(model, id) { + AKANTU_DEBUG_IN(); + + E1 = 0; + E2 = 0; + E3 = 0; + nu12 = 0; + nu13 = 0; + nu23 = 0; + G12 = 0; + G13 = 0; + G23 = 0; + + plane_stress = false; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialElasticOrthotropic::initMaterial() { + AKANTU_DEBUG_IN(); + Material::initMaterial(); + + S = new Real[9]; + Real Delta; + + Real nu21 = nu12 * E2 / E1; + Real nu31 = nu13 * E3 / E1; + Real nu32 = nu23 * E3 / E2; + + if(plane_stress) { + Delta = 1 - nu12 * nu21; + + S[0] = E1 / Delta; + S[1] = E2 / Delta; + S[2] = 0; + S[3] = E1 * nu21 / Delta; + S[4] = E2 * nu12 / Delta; + S[5] = 0; + S[6] = G12; + S[7] = 0; + S[8] = 0; + } + else { + Delta = 1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu13 * nu32; + + S[0] = E1 * (1 - nu23 * nu32) / Delta; + S[1] = E2 * (1 - nu31 * nu13) / Delta; + S[2] = E3 * (1 - nu12 * nu21) / Delta; + + S[3] = E2 * (nu12 + nu32 * nu13) / Delta; + S[4] = E1 * (nu31 + nu21 * nu32) / Delta; + S[5] = E3 * (nu23 + nu21 * nu13) / Delta; + + S[6] = G12; + S[7] = G13; + S[8] = G23; + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialElasticOrthotropic::computeStress(ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + Real F[3*3]; + Real sigma[3*3]; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; + memset(F, 0, 3 * 3 * sizeof(Real)); + + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + F[3*i + j] = strain_val[spatial_dimension * i + j]; + + // for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1; + + computeStress(F, sigma); + + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + stress_val[spatial_dimension*i + j] = sigma[3 * i + j]; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialElasticOrthotropic::computeTangentStiffness(const ElementType & el_type, + Vector & tangent_matrix, + GhostType ghost_type) { + + switch(spatial_dimension) { + case 1: { computeTangentStiffnessByDim<1>(el_type, tangent_matrix, ghost_type); break; } + case 2: { computeTangentStiffnessByDim<2>(el_type, tangent_matrix, ghost_type); break; } + case 3: { computeTangentStiffnessByDim<3>(el_type, tangent_matrix, ghost_type); break; } + } +} + +/* -------------------------------------------------------------------------- */ +bool MaterialElasticOrthotropic::setParam(const std::string & key, const std::string & value, + const ID & id) { + std::stringstream sstr(value); + if(key == "E1") { sstr >> E1; } + else if(key == "E2") { sstr >> E2; } + else if(key == "E3") { sstr >> E3; } + else if(key == "nu12") { sstr >> nu12; } + else if(key == "nu13") { sstr >> nu13; } + else if(key == "nu23") { sstr >> nu23; } + else if(key == "G12") { sstr >> G12; } + else if(key == "G13") { sstr >> G13; } + else if(key == "G23") { sstr >> G23; } + else if(key == "Plane_Stress") { sstr >> plane_stress; } + else { return Material::setParam(key, value, id); } + return true; +} + +/* -------------------------------------------------------------------------- */ +Real MaterialElasticOrthotropic::getPushWaveSpeed() { + return sqrt( (*std::max_element(S, S+3)) /rho); +} + +/* -------------------------------------------------------------------------- */ +Real MaterialElasticOrthotropic::getShearWaveSpeed() { + return sqrt( (*std::max_element(S+6, S+9)) /rho); +} + +/* -------------------------------------------------------------------------- */ +void MaterialElasticOrthotropic::printself(std::ostream & stream, int indent) const { + std::string space; + for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + + stream << space << "Material<_elastic<_orthotropic> [" << std::endl; + if(!plane_stress) + stream << space << " + Plane strain" << std::endl; + else + stream << space << " + Plane stress" << std::endl; + stream << space << " + density : " << rho << std::endl; + stream << space << " + Young's modulus (x) : " << E1 << std::endl; + stream << space << " + Young's modulus (y) : " << E2 << std::endl; + stream << space << " + Young's modulus (z) : " << E3 << std::endl; + stream << space << " + Poisson's ratio (xy) : " << nu12 << std::endl; + stream << space << " + Poisson's ratio (xz) : " << nu13 << std::endl; + stream << space << " + Poisson's ratio (yz) : " << nu23 << std::endl; + stream << space << " + Shear modulus (xy) : " << G12 << std::endl; + stream << space << " + Shear modulus (xz) : " << G13 << std::endl; + stream << space << " + Shear modulus (yz) : " << G23 << std::endl; + + Material::printself(stream, indent + 1); + stream << space << "]" << std::endl; +} +/* -------------------------------------------------------------------------- */ + +__END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh b/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh new file mode 100644 index 000000000..c237eafde --- /dev/null +++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.hh @@ -0,0 +1,181 @@ +/** + * @file material_elastic_orthotropic.hh + * @author Marco Vocialta + * @date Thu Apr 12 10:57:52 2012 + * + * @brief Orthotropic elastic material + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "material.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__ +#define __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__ + +__BEGIN_AKANTU__ + +/** + * Orthotropic elastic material + * + * parameters in the material files : + * - rho : density (default: 0) + * - E1 : Young's modulus along x (default: 0) + * - E2 : Young's modulus along y (default: 0) + * - E3 : Young's modulus along z (default: 0) + * - nu12 : Poisson's ratio along xy (default: 0) + * - nu13 : Poisson's ratio along xz (default: 0) + * - nu23 : Poisson's ratio along yz (default: 0) + * - G12 : Shear modulus along xy (default: 0) + * - G13 : Shear modulus along xz (default: 0) + * - G23 : Shear modulus along yz (default: 0) + * - Plane_Stress : if 0: plane strain, else: plane stress (default: 0) + */ +class MaterialElasticOrthotropic : public virtual Material { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +public: + + MaterialElasticOrthotropic(Model & model, const ID & id = ""); + + virtual ~MaterialElasticOrthotropic() {}; + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + + virtual void initMaterial(); + + virtual bool setParam(const std::string & key, const std::string & value, + const ID & id); + + /// constitutive law for all element of a type + virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); + + /// compute the tangent stiffness matrix for an element type + void computeTangentStiffness(const ElementType & el_type, + Vector & tangent_matrix, + GhostType ghost_type = _not_ghost); + + /// compute the p-wave speed in the material + virtual Real getPushWaveSpeed(); + + /// compute the s-wave speed in the material + virtual Real getShearWaveSpeed(); + + /// function to print the containt of the class + virtual void printself(std::ostream & stream, int indent = 0) const; + +protected: + /// constitutive law for a given quadrature point + inline void computeStress(Real * F, Real * sigma); + + // /// compute the tangent stiffness matrix for an element type + template + void computeTangentStiffnessByDim(akantu::ElementType, akantu::Vector& tangent_matrix, akantu::GhostType); + + // /// compute the tangent stiffness matrix for an element + template + void computeTangentStiffness(Real * tangent); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + /// get the stable time step + inline Real getStableTimeStep(Real h, const Element & element); + + AKANTU_GET_MACRO(E1, E1, Real); + AKANTU_GET_MACRO(E2, E2, Real); + AKANTU_GET_MACRO(E3, E3, Real); + AKANTU_GET_MACRO(Nu12, nu12, Real); + AKANTU_GET_MACRO(Nu13, nu13, Real); + AKANTU_GET_MACRO(Nu23, nu23, Real); + AKANTU_GET_MACRO(G12, G12, Real); + AKANTU_GET_MACRO(G13, G13, Real); + AKANTU_GET_MACRO(G23, G23, Real); + + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +protected: + + /// the x young modulus + Real E1; + + /// the y young modulus + Real E2; + + /// the z young modulus + Real E3; + + /// xy Poisson coefficient + Real nu12; + + /// xz Poisson coefficient + Real nu13; + + /// yz Poisson coefficient + Real nu23; + + /// xy shear modulus + Real G12; + + /// xz shear modulus + Real G13; + + /// yz shear modulus + Real G23; + + /// stiffness coefficients + Real * S; + + /// Plane stress or plane strain + bool plane_stress; + +}; + +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ + +#if defined (AKANTU_INCLUDE_INLINE_IMPL) +# include "material_elastic_orthotropic_inline_impl.cc" +#endif + +/* -------------------------------------------------------------------------- */ +/// standard output stream operator +inline std::ostream & operator <<(std::ostream & stream, const MaterialElasticOrthotropic & _this) +{ + _this.printself(stream); + return stream; +} + +__END_AKANTU__ + +#endif /* __AKANTU_MATERIAL_ELASTIC_ORTHOTROPIC_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic_inline_impl.cc b/src/model/solid_mechanics/materials/material_elastic_orthotropic_inline_impl.cc new file mode 100644 index 000000000..61169d63e --- /dev/null +++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic_inline_impl.cc @@ -0,0 +1,98 @@ +/** + * @file material_elastic_orthotropic_inline_impl.cc + * @author Marco Vocialta + * @date Thu Apr 12 13:40:42 2012 + * + * @brief Implementation of the inline functions of the orthotropic + * elastic material + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ + + +/* -------------------------------------------------------------------------- */ +inline void MaterialElasticOrthotropic::computeStress(Real * F, Real * sigma) { + /// \mathbf{\sigma} = \mathbf{S} \mathbf{F} + sigma[0] = S[0] * F[0] + S[3] * F[4] + S[4] * F[8]; + sigma[4] = S[3] * F[0] + S[1] * F[4] + S[5] * F[8]; + sigma[8] = S[4] * F[0] + S[5] * F[4] + S[2] * F[8]; + + sigma[1] = sigma[3] = S[6] * (F[1] + F[3]) / 2; + sigma[2] = sigma[6] = S[7] * (F[2] + F[6]) / 2; + sigma[5] = sigma[7] = S[8] * (F[5] + F[7]) / 2; +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialElasticOrthotropic::computeTangentStiffnessByDim(__attribute__((unused)) akantu::ElementType, + akantu::Vector& tangent_matrix, + __attribute__((unused)) akantu::GhostType) { + AKANTU_DEBUG_IN(); + + Real * tangent_val = tangent_matrix.values; + UInt offset_tangent = tangent_matrix.getNbComponent(); + UInt nb_quads = tangent_matrix.getSize(); + + if (nb_quads == 0) return; + + memset(tangent_val, 0, offset_tangent * nb_quads * sizeof(Real)); + for (UInt q = 0; q < nb_quads; ++q, tangent_val += offset_tangent) { + computeTangentStiffness(tangent_val); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialElasticOrthotropic::computeTangentStiffness(Real * tangent) { + + UInt n = (dim * (dim - 1) / 2 + dim); + + tangent[0 * n + 0] = S[0]; + + // test of dimension should by optimized out by the compiler due to the template + if(dim >= 2) { + tangent[1 * n + 1] = S[1]; + tangent[0 * n + 1] = S[3]; + tangent[1 * n + 0] = S[3]; + + tangent[(n - 1) * n + (n - 1)] = S[6]; + } + + if(dim == 3) { + tangent[2 * n + 2] = S[2]; + tangent[0 * n + 2] = S[4]; + tangent[1 * n + 2] = S[5]; + tangent[2 * n + 0] = S[4]; + tangent[2 * n + 1] = S[5]; + + tangent[3 * n + 3] = S[7]; + tangent[4 * n + 4] = S[8]; + } +} + +/* -------------------------------------------------------------------------- */ +inline Real MaterialElasticOrthotropic::getStableTimeStep(Real h, + __attribute__ ((unused)) const Element & element) { + return (h/getPushWaveSpeed()); +} diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc new file mode 100644 index 000000000..f35561ad7 --- /dev/null +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc @@ -0,0 +1,676 @@ +/** + * @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) { + 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; + MaterialCohesive * mat_cohesive; + + while ((mat_cohesive = + 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(std::string material_file, + AnalysisMethod method) { + + AKANTU_DEBUG_IN(); + + SolidMechanicsModelCohesive::initFull(material_file, method); + + MaterialCohesiveLinearExtrinsic * mat_cohesive + = dynamic_cast(materials[cohesive_index]); + AKANTU_DEBUG_ASSERT(mat_cohesive, "Choose the linear extrinsic cohesive constitutive law"); + + 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)); + + for (UInt f = 0; f < nb_facet; ++f) + sigma_lim(f) = sigma_c * (1 + std::rand()/(Real)RAND_MAX * rand); + + registerFEMObject("FacetsFEM", mesh_facets, spatial_dimension-1); + getFEM("FacetsFEM").initShapeFunctions(); + getFEM("FacetsFEM").computeNormalsOnControlPoints(); + + AKANTU_DEBUG_OUT(); + +} + +/* -------------------------------------------------------------------------- */ + +void SolidMechanicsModelCohesive::checkCohesiveStress() { + AKANTU_DEBUG_IN(); + + UInt nb_facet = mesh_facets.getNbElement(type_facet); + UInt nb_quad_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); + Vector facet_insertion; + Vector sigma_insertion; + + const Vector & normals = getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet); + + const Vector > & element_to_facet + = mesh_facets.getElementToSubelement(type_facet); + + for (UInt f = 0; f < nb_facet; ++f) { + /// check that facet isn't on the boundary + if (element_to_facet(f)(1).type != _not_defined) { + + Vector stress_on_quad(2*nb_quad_facet, spatial_dimension); + + for (UInt el = 0; el < 2; ++el) { + UInt elem_nb = element_to_facet(f)(el).element; + ElementType type_elem = element_to_facet(f)(el).type; + // UInt nb_quad_elem = getFEM().getNbQuadraturePoints(type_elem); + const Vector & stress = getMaterial(0).getStress(type_elem); + + /// linear element case + types::Matrix stress_local(spatial_dimension, spatial_dimension); + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + stress_local(i, j) = stress(elem_nb, i*spatial_dimension + j); + + types::RVector normal(spatial_dimension); + for (UInt i = 0; i < spatial_dimension; ++i) + normal(i) = normals(f, i); + + types::RVector stress_norm(spatial_dimension); + stress_norm.mul(stress_local, normal); + + for (UInt dim = 0; dim < spatial_dimension; ++dim) + stress_on_quad(el, dim) = fabs(stress_norm(dim)); + } + + Real sigma_max = *std::max_element(stress_on_quad.storage(), + stress_on_quad.storage() + + stress_on_quad.getSize() * + stress_on_quad.getNbComponent()); + + if (sigma_lim(f) < sigma_max) { + facet_insertion.push_back(f); + sigma_insertion.push_back(sigma_max); + } + } + } + + insertCohesiveElements(facet_insertion, true); + + /// resize cohesive material vectors + MaterialCohesiveLinearExtrinsic * mat_cohesive + = dynamic_cast(materials[cohesive_index]); + AKANTU_DEBUG_ASSERT(mat_cohesive, "Initialize the material with initExtrinsic"); + mat_cohesive->resizeCohesiveVectors(sigma_insertion); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ + +void SolidMechanicsModelCohesive::insertCohesiveElements(const Vector & facet_insertion, + const bool random_sigma) { + + AKANTU_DEBUG_IN(); + + Vector doubled_nodes(0, 2); + Vector doubled_facets(0, 2); + + /// update mesh + for (UInt f = 0; f < facet_insertion.getSize(); ++f) { + Element facet; + facet.element = facet_insertion(f); + facet.type = type_facet; + doubleFacet(facet, doubled_nodes, doubled_facets); + } + + /// double middle nodes if it's the case + if (type_facet == _segment_3) { + doubleMiddleNode(doubled_nodes, doubled_facets); + } + + /// update nodal values + updateDoubledNodes(doubled_nodes); + + /// 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); + + 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 shape functions + getFEM("CohesiveFEM").initShapeFunctions(_not_ghost); + + /// resize cohesive material vectors + if (!random_sigma) { + MaterialCohesive * mat_cohesive + = dynamic_cast(materials[cohesive_index]); + AKANTU_DEBUG_ASSERT(mat_cohesive, "No cohesive materials in the materials vector"); + mat_cohesive->resizeCohesiveVectors(); + } + + 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 + element_to_facet.push_back(element_to_facet(f_index)); + + /// 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) { + 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; + current_facet.element = f_global; + current_facet.type = type_facet; + 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::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 new file mode 100644 index 000000000..e06151d52 --- /dev/null +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh @@ -0,0 +1,161 @@ +/** + * @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, + const bool random_sigma = false); + + /// initialize completely the model + void initFull(std::string material_file, + AnalysisMethod method = _explicit_dynamic); + + /// initialize completely the model for extrinsic elements + void initExtrinsic(std::string material_file, + AnalysisMethod method = _explicit_dynamic); + + /// initialize the model + void initModel(); + + /// initialize cohesive material + void initCohesiveMaterial(); + +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 facet type + AKANTU_GET_MACRO(FacetType, type_facet, ElementType); + + /// get the facet mesh + AKANTU_GET_MACRO(MeshFacets, mesh_facets, 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 &); + + /* ------------------------------------------------------------------------ */ + /* 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; + +}; + + +/* -------------------------------------------------------------------------- */ +/* 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_cohesive/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_cohesive/CMakeLists.txt new file mode 100644 index 000000000..56705247f --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/CMakeLists.txt @@ -0,0 +1,48 @@ +#=============================================================================== +# @file CMakeLists.txt +# @author Marco Vocialta +# @date Fri Feb 24 14:30:59 2012 +# +# @section LICENSE +# +# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . +# +# @section DESCRIPTION +# +#=============================================================================== + +add_akantu_test(test_cohesive_buildfacets "test_cohesive_buildfacets") +add_akantu_test(test_cohesive_intrinsic "test_cohesive_intrinsic") +add_akantu_test(test_cohesive_extrinsic "test_cohesive_extrinsic") + +#=============================================================================== + +register_test(test_cohesive test_cohesive.cc) + +#=============================================================================== +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/material.dat + ${CMAKE_CURRENT_BINARY_DIR}/material.dat + COPYONLY + ) +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/mesh.msh + ${CMAKE_CURRENT_BINARY_DIR}/mesh.msh + COPYONLY + ) +file(COPY run DESTINATION .) +file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/material.dat new file mode 100644 index 000000000..18c7bd799 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/material.dat @@ -0,0 +1,15 @@ +material elastic [ + name = steel + rho = 1e7 # density + E = 1e7 # young's modulus + nu = 0.3 # poisson's ratio +] + +material cohesive_linear_extrinsic [ + name = cohesive + sigma_c = 1e5 + beta = 0 + G_cI = 1e5 + G_cII = 1e5 + rand = 20 +] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/mesh.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh.geo new file mode 100644 index 000000000..aa2462049 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh.geo @@ -0,0 +1,18 @@ +h =0.2; + +Point(1) = {0, 0, 0, h}; +Point(2) = {0, 1, 0, h}; +Point(3) = {-1, 1, 0, h}; +Point(4) = {-1, 0, 0, h}; +Point(6) = {1, 0, 0, h}; +Point(7) = {1, 1, 0, h}; +Line(1) = {2, 3}; +Line(2) = {3, 4}; +Line(3) = {4, 1}; +Line(4) = {1, 6}; +Line(5) = {6, 7}; +Line(6) = {7, 2}; + +Line Loop(7) = {1, 2, 3, 4, 5, 6}; +Plane Surface(8) = {7}; +Physical Surface(0) = {8}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/mesh.msh b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh.msh new file mode 100644 index 000000000..8cdab28b1 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh.msh @@ -0,0 +1,142 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +53 +1 0 0 0 +2 0 1 0 +3 -1 1 0 +4 -1 0 0 +5 1 0 0 +6 1 1 0 +7 -0.2499999999993998 1 0 +8 -0.4999999999990217 1 0 +9 -0.7499999999993374 1 0 +10 -1 0.7500000000003469 0 +11 -1 0.5000000000013878 0 +12 -1 0.2500000000010407 0 +13 -0.7500000000003469 0 0 +14 -0.5000000000013878 0 0 +15 -0.2500000000010407 0 0 +16 0.2499999999993998 0 0 +17 0.4999999999990217 0 0 +18 0.7499999999993374 0 0 +19 1 0.2499999999993998 0 +20 1 0.4999999999990217 0 +21 1 0.7499999999993374 0 +22 0.7500000000003469 1 0 +23 0.5000000000013878 1 0 +24 0.2500000000010407 1 0 +25 -0.109885838682046 0.4966964540720109 0 +26 0.3855511280698393 0.5000000000003061 0 +27 -0.5501220340404019 0.4080683147660409 0 +28 0.6510655547138771 0.3176045927728005 0 +29 0.6510655547140901 0.6823954072276218 0 +30 0.1125309742557222 0.2618727640254615 0 +31 0.1089015214619986 0.7358134564467858 0 +32 -0.3574894723948559 0.7281878845348634 0 +33 -0.7266746693214478 0.7019059774796227 0 +34 -0.3341676355493891 0.2692973568512923 0 +35 0.4025972431137625 0.7356032509683735 0 +36 0.4028996975127416 0.2642039340714447 0 +37 -0.1219506414480373 0.7838786635824847 0 +38 -0.7450297577945056 0.252179537794927 0 +39 -0.1074687777815134 0.2098576111188913 0 +40 -0.7889741573686897 0.4960941618476821 0 +41 -0.5596752946346957 0.1966354880432951 0 +42 -0.371703150536428 0.4979169684315345 0 +43 0.8062022896298995 0.4999999999999543 0 +44 0.205579736632924 0.4993544288428088 0 +45 0.8080003462072091 0.8032599668670146 0 +46 0.8080003462067507 0.1967400331328585 0 +47 0.5610725594392284 0.4999785761156746 0 +48 -0.591984853628037 0.8123632048735649 0 +49 0.6232620797906236 0.1514034669779444 0 +50 0.6232217525385286 0.8485708243608364 0 +51 -0.8206602848780158 0.8467415727987578 0 +52 -0.5564796491946886 0.611350844142008 0 +53 -0.8721378961405318 0.1273969864581581 0 +$EndNodes +$Elements +80 +1 2 2 0 8 1 16 30 +2 2 2 0 8 2 31 24 +3 2 2 0 8 7 8 32 +4 2 2 0 8 8 48 32 +5 2 2 0 8 14 15 34 +6 2 2 0 8 10 40 33 +7 2 2 0 8 14 34 41 +8 2 2 0 8 10 11 40 +9 2 2 0 8 2 37 31 +10 2 2 0 8 19 43 46 +11 2 2 0 8 21 45 43 +12 2 2 0 8 16 36 30 +13 2 2 0 8 24 31 35 +14 2 2 0 8 25 31 37 +15 2 2 0 8 25 30 44 +16 2 2 0 8 25 44 31 +17 2 2 0 8 7 32 37 +18 2 2 0 8 25 37 32 +19 2 2 0 8 29 35 47 +20 2 2 0 8 28 47 36 +21 2 2 0 8 25 39 30 +22 2 2 0 8 26 35 44 +23 2 2 0 8 31 44 35 +24 2 2 0 8 26 44 36 +25 2 2 0 8 30 36 44 +26 2 2 0 8 12 38 40 +27 2 2 0 8 11 12 40 +28 2 2 0 8 1 30 39 +29 2 2 0 8 26 47 35 +30 2 2 0 8 26 36 47 +31 2 2 0 8 19 20 43 +32 2 2 0 8 20 21 43 +33 2 2 0 8 25 32 42 +34 2 2 0 8 28 46 43 +35 2 2 0 8 29 43 45 +36 2 2 0 8 10 33 51 +37 2 2 0 8 13 41 38 +38 2 2 0 8 25 34 39 +39 2 2 0 8 13 14 41 +40 2 2 0 8 15 39 34 +41 2 2 0 8 9 51 48 +42 2 2 0 8 16 17 36 +43 2 2 0 8 23 24 35 +44 2 2 0 8 29 50 35 +45 2 2 0 8 28 36 49 +46 2 2 0 8 33 48 51 +47 2 2 0 8 2 7 37 +48 2 2 0 8 25 42 34 +49 2 2 0 8 27 41 34 +50 2 2 0 8 1 39 15 +51 2 2 0 8 27 34 42 +52 2 2 0 8 8 9 48 +53 2 2 0 8 23 35 50 +54 2 2 0 8 17 49 36 +55 2 2 0 8 5 46 18 +56 2 2 0 8 6 22 45 +57 2 2 0 8 5 19 46 +58 2 2 0 8 6 45 21 +59 2 2 0 8 27 40 38 +60 2 2 0 8 27 38 41 +61 2 2 0 8 27 52 40 +62 2 2 0 8 33 40 52 +63 2 2 0 8 3 51 9 +64 2 2 0 8 3 10 51 +65 2 2 0 8 13 38 53 +66 2 2 0 8 12 53 38 +67 2 2 0 8 17 18 49 +68 2 2 0 8 22 23 50 +69 2 2 0 8 4 53 12 +70 2 2 0 8 4 13 53 +71 2 2 0 8 27 42 52 +72 2 2 0 8 42 32 52 +73 2 2 0 8 28 43 47 +74 2 2 0 8 29 47 43 +75 2 2 0 8 28 49 46 +76 2 2 0 8 18 46 49 +77 2 2 0 8 22 50 45 +78 2 2 0 8 29 45 50 +79 2 2 0 8 32 48 52 +80 2 2 0 8 33 52 48 +$EndElements diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_3d.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_3d.geo new file mode 100644 index 000000000..b46bc16df --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_3d.geo @@ -0,0 +1,35 @@ +Point(1) = {0, 0, 0, 1.0}; +Point(2) = {1, 0, 0, 1.0}; +Point(3) = {0, 1, 0, 1.0}; +Point(4) = {0, -1, 0, 1.0}; +Point(5) = {-1, 0, 0, 1.0}; +Point(7) = {0.5, 0.5, 0, 1.0}; +Point(8) = {0.5, -0.5, 0, 1.0}; +Point(9) = {-0.5, 0.5, 0, 1.0}; +Point(10) = {-0.5, -0.5, 0, 1.0}; +Point(11) = {0, 0, 1, 1.0}; +Line(1) = {4, 8}; +Line(2) = {8, 2}; +Line(3) = {2, 7}; +Line(4) = {7, 3}; +Line(5) = {3, 9}; +Line(6) = {9, 5}; +Line(7) = {5, 10}; +Line(8) = {10, 4}; +Line(9) = {4, 11}; +Line(10) = {11, 5}; +Line(11) = {11, 3}; +Line(12) = {11, 2}; +Line Loop(13) = {2, 3, 4, 5, 6, 7, 8, 1}; +Plane Surface(14) = {13}; +Line Loop(15) = {10, 7, 8, 9}; +Plane Surface(16) = {15}; +Line Loop(17) = {9, 12, -2, -1}; +Plane Surface(18) = {17}; +Line Loop(19) = {12, 3, 4, -11}; +Plane Surface(20) = {19}; +Line Loop(21) = {11, 5, 6, -10}; +Plane Surface(22) = {21}; +Surface Loop(23) = {18, 16, 22, 20, 14}; +Volume(24) = {23}; +Physical Volume(25) = {24}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old.geo new file mode 100644 index 000000000..cd4bb3434 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old.geo @@ -0,0 +1,25 @@ +h = 1.0; + +Point(1) = {0, 0, 0, h}; +Point(2) = {0, 1, 0, h}; +Point(3) = {-1, 1, 0, h}; +Point(4) = {-1, 0, 0, h}; +Point(5) = {0, 0, 0, h}; +Point(6) = {1, 0, 0, h}; +Point(7) = {1, 1, 0, h}; +Point(8) = {0, 1, 0, h}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line(5) = {5, 6}; +Line(6) = {6, 7}; +Line(7) = {7, 8}; +Line(8) = {8, 5}; + +Line Loop(9) = {5, 6, 7, 8}; +Plane Surface(10) = {9}; +Line Loop(11) = {1, 2, 3, 4}; +Plane Surface(12) = {11}; +Physical Surface(0) = {12}; +Physical Surface(1) = {10}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old.msh b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old.msh new file mode 100644 index 000000000..6d3952bd7 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old.msh @@ -0,0 +1,27 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +10 +1 0 0 0 +2 0 1 0 +3 -1 1 0 +4 -1 0 0 +5 0 0 0 +6 1 0 0 +7 1 1 0 +8 0 1 0 +9 0.5 0.5 0 +10 -0.5 0.5 0 +$EndNodes +$Elements +8 +1 2 2 14 10 5 9 8 +2 2 2 14 10 5 6 9 +3 2 2 14 10 6 7 9 +4 2 2 14 10 7 8 9 +5 2 2 13 12 3 4 10 +6 2 2 13 12 2 3 10 +7 2 2 13 12 2 10 1 +8 2 2 13 12 4 1 10 +$EndElements diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old3d.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old3d.geo new file mode 100644 index 000000000..56dfe7ac3 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/mesh_old3d.geo @@ -0,0 +1,24 @@ +lc = 1.0; +Point(1) = {0, 0, 0, lc}; +Point(2) = {1, 0, 0, lc}; +Point(3) = {0, 1, 0, lc}; +Point(4) = {-1, 0, 0, lc}; +Point(5) = {0, 0, 1, lc}; +Line(1) = {1, 2}; +Line(3) = {2, 3}; +Line(4) = {3, 4}; +Line(5) = {4, 1}; +Line(6) = {4, 5}; +Line(7) = {5, 2}; +Line(8) = {5, 3}; +Line Loop(9) = {1, 3, 4, 5}; +Plane Surface(10) = {9}; +Line Loop(11) = {4, 6, 8}; +Plane Surface(12) = {11}; +Line Loop(13) = {3, -8, 7}; +Plane Surface(14) = {13}; +Line Loop(15) = {6, 7, -1, -5}; +Plane Surface(16) = {15}; +Surface Loop(17) = {16, 12, 10, 14}; +Volume(18) = {17}; +Physical Volume(19) = {18}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/run b/test/test_model/test_solid_mechanics_model/test_cohesive/run new file mode 100755 index 000000000..420dc952e --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/run @@ -0,0 +1,7 @@ +#! /bin/bash + +rm -fr paraview +mkdir paraview + +./test_cohesive + diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive.cc new file mode 100644 index 000000000..289d86bde --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive.cc @@ -0,0 +1,202 @@ +/** + * @file test_cohesive.cc + * @author Marco Vocialta + * @date Fri Feb 24 14:32:31 2012 + * + * @brief Test 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 "aka_common.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "material.hh" +#include "io_helper.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + initialize(argc, argv); + + debug::setDebugLevel(dblWarning); + + const UInt spatial_dimension = 2; + const UInt max_steps = 2500; + + const ElementType type = _triangle_3; + + Mesh mesh(spatial_dimension); + MeshIOMSH mesh_io; + mesh_io.read("mesh.msh", mesh); + + SolidMechanicsModelCohesive model(mesh); + + /// model initialization + model.initExtrinsic("material.dat"); + Real time_step = model.getStableTimeStep()*0.8; + model.setTimeStep(time_step); + std::cout << "Time step: " << time_step << std::endl; + + model.assembleMassLumped(); + + const ElementType type_facet = mesh.getFacetElementType(type); + + const Mesh & mesh_facets = model.getMeshFacets(); + + Vector & boundary = model.getBoundary(); + // Vector & velocity = model.getVelocity(); + Vector & displacement = model.getDisplacement(); + const Vector & position = mesh.getNodes(); + // Vector & acceleration = model.getAcceleration(); + // Vector & increment_acceleration = model.getIncrementAcceleration(); + // Vector & force = model.getResidual(); + // const Vector & strain = model.getMaterial(0).getStrain(type); + // const Vector & stress = model.getMaterial(0).getStress(type); + + UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); + + /// boundary conditions + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 0) == -1 || position(n, 0) == 1) { + boundary(n, 0) = boundary(n, 1) = true; + } + } + + + iohelper::ElemType paraview_type = iohelper::TRIANGLE1; + UInt nb_element = model.getFEM().getMesh().getNbElement(type); + + model.updateResidual(); + + /// initialize the paraview output + iohelper::DumperParaview dumper; + // dumper.SetMode(iohelper::TEXT); + dumper.SetPoints(model.getFEM().getMesh().getNodes().values, + spatial_dimension, mesh.getNbNodes(), "explicit"); + dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values, + paraview_type, nb_element, iohelper::C_MODE); + dumper.AddNodeDataField(model.getDisplacement().values, + spatial_dimension, "displacements"); + dumper.AddNodeDataField(model.getVelocity().values, + spatial_dimension, "velocity"); + dumper.AddNodeDataField(model.getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddNodeDataField(model.getForce().values, + spatial_dimension, "applied_force"); + dumper.AddNodeDataField(model.getResidual().values, + spatial_dimension, "forces"); + dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + spatial_dimension*spatial_dimension, "stress"); + dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + spatial_dimension*spatial_dimension, "strain"); + dumper.SetEmbeddedValue("displacements", 1); + dumper.SetEmbeddedValue("applied_force", 1); + dumper.SetEmbeddedValue("forces", 1); + dumper.SetPrefix("paraview/"); + dumper.Init(); + dumper.Dump(); + + + /// update displacement + + Real increment = 0.0001; + + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 0) == 1) + displacement(n, 0) += increment; + } + + + std::ofstream edis("edis.txt"); + std::ofstream erev("erev.txt"); + + /// main loop + for (UInt s = 1; s <= max_steps; ++s) { + + model.checkCohesiveStress(); + + model.explicitPred(); + model.updateResidual(); + model.updateAcceleration(); + model.explicitCorr(); + + if(s % 1 == 0) { + + dumper.SetPoints(model.getFEM().getMesh().getNodes().values, + spatial_dimension, mesh.getNbNodes(), "explicit"); + dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values, + paraview_type, nb_element, iohelper::C_MODE); + dumper.AddNodeDataField(model.getDisplacement().values, + spatial_dimension, "displacements"); + dumper.AddNodeDataField(model.getVelocity().values, + spatial_dimension, "velocity"); + dumper.AddNodeDataField(model.getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddNodeDataField(model.getForce().values, + spatial_dimension, "applied_force"); + dumper.AddNodeDataField(model.getResidual().values, + spatial_dimension, "forces"); + dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + spatial_dimension*spatial_dimension, "strain"); + dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + spatial_dimension*spatial_dimension, "stress"); + dumper.SetEmbeddedValue("displacements", 1); + dumper.SetEmbeddedValue("applied_force", 1); + dumper.SetEmbeddedValue("forces", 1); + dumper.SetPrefix("paraview/"); + dumper.Init(); + dumper.Dump(); + std::cout << "passing step " << s << "/" << max_steps << std::endl; + } + + /// update displacement + for (UInt n = 0; n < nb_nodes; ++n) { + if (position(n, 0) == 1) + displacement(n, 0) += increment; + } + + Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); + + edis << s << " " + << Ed << std::endl; + + erev << s << " " + << Er << std::endl; + + } + + edis.close(); + erev.close(); + + finalize(); + return EXIT_SUCCESS; +} diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/CMakeLists.txt new file mode 100644 index 000000000..21b990328 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/CMakeLists.txt @@ -0,0 +1,33 @@ +#=============================================================================== +# @file CMakeLists.txt +# @author Marco Vocialta +# @date Wed Apr 18 17:16:30 2012 +# +# @section LICENSE +# +# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . +# +# @section DESCRIPTION +# +#=============================================================================== +register_test(test_cohesive_buildfacets test_cohesive_buildfacets.cc) +file(COPY mesh.msh DESTINATION .) + +#add_mesh(test_cohesive_buildfacets_mesh mesh.geo 3 2) + +#add_dependencies(test_cohesive_buildfacets test_cohesive_buildfacets_mesh) + diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/mesh.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/mesh.geo new file mode 100644 index 000000000..b46bc16df --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/mesh.geo @@ -0,0 +1,35 @@ +Point(1) = {0, 0, 0, 1.0}; +Point(2) = {1, 0, 0, 1.0}; +Point(3) = {0, 1, 0, 1.0}; +Point(4) = {0, -1, 0, 1.0}; +Point(5) = {-1, 0, 0, 1.0}; +Point(7) = {0.5, 0.5, 0, 1.0}; +Point(8) = {0.5, -0.5, 0, 1.0}; +Point(9) = {-0.5, 0.5, 0, 1.0}; +Point(10) = {-0.5, -0.5, 0, 1.0}; +Point(11) = {0, 0, 1, 1.0}; +Line(1) = {4, 8}; +Line(2) = {8, 2}; +Line(3) = {2, 7}; +Line(4) = {7, 3}; +Line(5) = {3, 9}; +Line(6) = {9, 5}; +Line(7) = {5, 10}; +Line(8) = {10, 4}; +Line(9) = {4, 11}; +Line(10) = {11, 5}; +Line(11) = {11, 3}; +Line(12) = {11, 2}; +Line Loop(13) = {2, 3, 4, 5, 6, 7, 8, 1}; +Plane Surface(14) = {13}; +Line Loop(15) = {10, 7, 8, 9}; +Plane Surface(16) = {15}; +Line Loop(17) = {9, 12, -2, -1}; +Plane Surface(18) = {17}; +Line Loop(19) = {12, 3, 4, -11}; +Plane Surface(20) = {19}; +Line Loop(21) = {11, 5, 6, -10}; +Plane Surface(22) = {21}; +Surface Loop(23) = {18, 16, 22, 20, 14}; +Volume(24) = {23}; +Physical Volume(25) = {24}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/mesh.msh b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/mesh.msh new file mode 100644 index 000000000..8dfc3bebf --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/mesh.msh @@ -0,0 +1,80 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +55 +1 1 0 0 +2 0 1 0 +3 0 -1 0 +4 -1 0 0 +5 0.5 0.5 0 +6 0.5 -0.5 0 +7 -0.5 0.5 0 +8 -0.5 -0.5 0 +9 0 0 1 +10 0.2499999999993232 -0.7500000000006768 0 +11 0.7499999999993227 -0.2500000000006773 0 +12 0.7500000000006768 0.2499999999993232 0 +13 0.2500000000006773 0.7499999999993227 0 +14 -0.2499999999993232 0.7500000000006768 0 +15 -0.7499999999993227 0.2500000000006773 0 +16 -0.7500000000006768 -0.2499999999993232 0 +17 -0.2500000000006773 -0.7499999999993227 0 +18 0 -0.5000000000013294 0.4999999999986707 +19 -0.4999999999986707 0 0.5000000000013294 +20 0 0.4999999999986707 0.5000000000013294 +21 0.4999999999986707 0 0.5000000000013294 +22 0 0 0 +23 0 -0.5 0 +24 -0.5 0 0 +25 0 0.5 0 +26 0.5 0 0 +27 -0.5 -0.25 0 +28 -0.75 0 0 +29 -0.25 0.5 0 +30 0 0.75 0 +31 0 -0.75 0 +32 -0.25 -0.5 0 +33 0.25 -0.5 0 +34 -0.5 0.25 0 +35 0.25 0.5 0 +36 0.75 0 0 +37 0.5 -0.25 0 +38 0.5 0.25 0 +39 -0.25 0 0 +40 -0.25 0.25 0 +41 0.25 -0.25 0 +42 0 -0.25 0 +43 -0.25 -0.25 0 +44 0 0.25 0 +45 0.25 0.25 0 +46 0.25 0 0 +47 -0.2500000000000001 -0.2499999999999999 0.5 +48 0.2499999999999999 -0.2500000000000002 0.5000000000000002 +49 0.2500000000000001 0.2499999999999999 0.5 +50 -0.2500000000000001 0.2499999999999999 0.5 +51 0 0 0.5 +52 0.25 0 0.5 +53 0 0.25 0.5 +54 0 -0.25 0.5 +55 -0.25 0 0.5 +$EndNodes +$Elements +16 +1 11 2 25 24 22 5 9 26 45 49 51 46 52 38 +2 11 2 25 24 22 9 25 7 51 53 44 40 29 50 +3 11 2 25 24 6 9 1 26 48 21 11 37 36 52 +4 11 2 25 24 9 23 3 8 54 31 18 47 17 32 +5 11 2 25 24 9 23 22 6 54 42 51 48 41 33 +6 11 2 25 24 23 9 3 6 54 18 31 33 10 48 +7 11 2 25 24 9 6 22 26 48 41 51 52 46 37 +8 11 2 25 24 24 9 8 22 55 47 27 39 43 51 +9 11 2 25 24 9 23 8 22 54 32 47 51 43 42 +10 11 2 25 24 9 24 8 4 55 27 47 19 16 28 +11 11 2 25 24 2 9 25 5 20 53 30 13 35 49 +12 11 2 25 24 22 5 25 9 45 35 44 51 53 49 +13 11 2 25 24 9 1 26 5 21 36 52 49 38 12 +14 11 2 25 24 24 9 22 7 55 51 39 34 40 50 +15 11 2 25 24 24 9 7 4 55 50 34 28 15 19 +16 11 2 25 24 9 2 25 7 20 30 53 50 29 14 +$EndElements diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/test_cohesive_buildfacets.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/test_cohesive_buildfacets.cc new file mode 100644 index 000000000..8f5864c3d --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfacets/test_cohesive_buildfacets.cc @@ -0,0 +1,253 @@ +/** + * @file test_cohesive.cc + * @author Marco Vocialta + * @date Fri Feb 24 14:32:31 2012 + * + * @brief Test 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 "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; + +int main(int argc, char *argv[]) { + initialize(argc, argv); + + // debug::setDebugLevel(dblDump); + + const UInt spatial_dimension = 3; + const ElementType type = _tetrahedron_10; + + Mesh mesh(spatial_dimension); + MeshIOMSH mesh_io; + mesh_io.read("mesh.msh", mesh); + Mesh mesh_facets(spatial_dimension, const_cast &>(mesh.getNodes()), "mesh_facets", 1); + + MeshUtils::buildAllFacets(mesh, mesh_facets); + + // std::cout << mesh << std::endl; + + const ElementType type_facet = mesh.getFacetElementType(type); + const ElementType type_subfacet = mesh.getFacetElementType(type_facet); + const ElementType type_subsubfacet = mesh.getFacetElementType(type_subfacet); + + /* ------------------------------------------------------------------------ */ + /* Element to Subelement testing */ + /* ------------------------------------------------------------------------ */ + + const Vector > & el_to_subel3 = mesh_facets.getElementToSubelement(type_facet); + const Vector > & el_to_subel2 = mesh_facets.getElementToSubelement(type_subfacet); + const Vector > & el_to_subel1 = mesh_facets.getElementToSubelement(type_subsubfacet); + + /// build vectors for comparison + Vector tetrahedron(2); + tetrahedron(0).type = type; + tetrahedron(0).element = 1; + + tetrahedron(1).type = type; + tetrahedron(1).element = 11; + + Vector triangle(8); + triangle(0).type = type_facet; + triangle(0).element = 2; + + triangle(1).type = type_facet; + triangle(1).element = 0; + + triangle(2).type = type_facet; + triangle(2).element = 4; + + triangle(3).type = type_facet; + triangle(3).element = 7; + + triangle(4).type = type_facet; + triangle(4).element = 26; + + triangle(5).type = type_facet; + triangle(5).element = 24; + + triangle(6).type = type_facet; + triangle(6).element = 16; + + triangle(7).type = type_facet; + triangle(7).element = 18; + + Vector segment(13); + segment(0).type = type_subfacet; + segment(0).element = 0; + + segment(1).type = type_subfacet; + segment(1).element = 1; + + segment(2).type = type_subfacet; + segment(2).element = 3; + + segment(3).type = type_subfacet; + segment(3).element = 7; + + segment(4).type = type_subfacet; + segment(4).element = 9; + + segment(5).type = type_subfacet; + segment(5).element = 12; + + segment(6).type = type_subfacet; + segment(6).element = 13; + + segment(7).type = type_subfacet; + segment(7).element = 16; + + segment(8).type = type_subfacet; + segment(8).element = 18; + + segment(9).type = type_subfacet; + segment(9).element = 21; + + segment(10).type = type_subfacet; + segment(10).element = 27; + + segment(11).type = type_subfacet; + segment(11).element = 32; + + segment(12).type = type_subfacet; + segment(12).element = 34; + + /// comparison + + for (UInt i = 0; i < tetrahedron.getSize(); ++i) { + if (tetrahedron(i).type != el_to_subel3(4)(i).type || + tetrahedron(i).element != el_to_subel3(4)(i).element) { + std::cout << tetrahedron(i).element << " " << el_to_subel3(4)(i).element << std::endl; + std::cout << "The two tetrahedrons connected to triangle 4 are wrong" + << std::endl; + return EXIT_FAILURE; + } + } + + for (UInt i = 0; i < triangle.getSize(); ++i) { + if (triangle(i).type != el_to_subel2(0)(i).type || + triangle(i).element != el_to_subel2(0)(i).element) { + std::cout << "The triangles connected to segment 0 are wrong" + << std::endl; + return EXIT_FAILURE; + } + } + + for (UInt i = 0; i < segment.getSize(); ++i) { + if (segment(i).type != el_to_subel1(1)(i).type || + segment(i).element != el_to_subel1(1)(i).element) { + std::cout << "The segments connected to point 1 are wrong" + << std::endl; + return EXIT_FAILURE; + } + } + + + /* ------------------------------------------------------------------------ */ + /* Subelement to Element testing */ + /* ------------------------------------------------------------------------ */ + + const Vector & subel_to_el3 = mesh_facets.getSubelementToElement(type); + const Vector & subel_to_el2 = mesh_facets.getSubelementToElement(type_facet); + const Vector & subel_to_el1 = mesh_facets.getSubelementToElement(type_subfacet); + + /// build vectors for comparison + Vector triangle2(mesh.getNbFacetsPerElement(type)); + triangle2(0).type = type_facet; + triangle2(0).element = 4; + + triangle2(1).type = type_facet; + triangle2(1).element = 5; + + triangle2(2).type = type_facet; + triangle2(2).element = 6; + + triangle2(3).type = type_facet; + triangle2(3).element = 7; + + Vector segment2(3); + segment2(0).type = type_subfacet; + segment2(0).element = 1; + + segment2(1).type = type_subfacet; + segment2(1).element = 3; + + segment2(2).type = type_subfacet; + segment2(2).element = 4; + + Vector point(2); + point(0).type = mesh.getFacetElementType(type_subfacet); + point(0).element = 1; + + point(1).type = mesh.getFacetElementType(type_subfacet); + point(1).element = 2; + + /// comparison + + for (UInt i = 0; i < triangle2.getSize(); ++i) { + if (triangle2(i).type != subel_to_el3(1, i).type || + triangle2(i).element != subel_to_el3(1, i).element) { + std::cout << "The triangles connected to tetrahedron 1 are wrong" + << std::endl; + return EXIT_FAILURE; + } + } + + for (UInt i = 0; i < segment2.getSize(); ++i) { + if (segment2(i).type != subel_to_el2(1, i).type || + segment2(i).element != subel_to_el2(1, i).element) { + std::cout << "The segments connected to triangle 1 are wrong" + << std::endl; + return EXIT_FAILURE; + } + } + + for (UInt i = 0; i < point.getSize(); ++i) { + if (point(i).type != subel_to_el1(1, i).type || + point(i).element != subel_to_el1(1, i).element) { + std::cout << "The points connected to segment 1 are wrong" + << std::endl; + return EXIT_FAILURE; + } + } + + + finalize(); + + std::cout << "OK: test_cohesive_buildfacets was passed!" << std::endl; + + return EXIT_SUCCESS; +} diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/CMakeLists.txt new file mode 100644 index 000000000..b19f49b07 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/CMakeLists.txt @@ -0,0 +1,39 @@ +#=============================================================================== +# @file CMakeLists.txt +# @author Marco Vocialta +# @date Thu Apr 26 15:14:00 2012 +# +# @section LICENSE +# +# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . +# +# @section DESCRIPTION +# +#=============================================================================== + +register_test(test_cohesive_extrinsic test_cohesive_extrinsic.cc) + +add_mesh(test_cohesive_extrinsic_mesh mesh.geo 2 1) +add_dependencies(test_cohesive_extrinsic test_cohesive_extrinsic_mesh) + +#=============================================================================== +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/material.dat + ${CMAKE_CURRENT_BINARY_DIR}/material.dat + COPYONLY + ) +#file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat new file mode 100644 index 000000000..9d6dfc631 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat @@ -0,0 +1,15 @@ +material elastic [ + name = steel + rho = 1e3 # density + E = 1e3 # young's modulus + nu = 0.001 # poisson's ratio +] + +material cohesive_linear_extrinsic [ + name = cohesive + sigma_c = 1 + beta = 0 + G_cI = 100 + G_cII = 100 + rand = 0 +] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo new file mode 100644 index 000000000..4526cbf34 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo @@ -0,0 +1,22 @@ +h = 1.0; + +Point(1) = {0, 0, 0, h}; +Point(2) = {0, 1, 0, h}; +Point(3) = {-1, 1, 0, h}; +Point(4) = {-1, 0, 0, h}; +Point(6) = {1, 0, 0, h}; +Point(7) = {1, 1, 0, h}; +Point(8) = {1, -1, 0, h}; +Point(9) = {0, -1, 0, h}; +Point(10) = {-1, -1, 0, h}; +Line(1) = {10, 9}; +Line(2) = {9, 8}; +Line(3) = {8, 6}; +Line(4) = {6, 7}; +Line(5) = {7, 2}; +Line(6) = {2, 3}; +Line(7) = {3, 4}; +Line(8) = {4, 10}; +Line Loop(9) = {3, 4, 5, 6, 7, 8, 1, 2}; +Plane Surface(10) = {9}; +Physical Surface(11) = {10}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc new file mode 100644 index 000000000..84a6afc75 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc @@ -0,0 +1,230 @@ +/** + * @file test_cohesive.cc + * @author Marco Vocialta + * @date Fri Feb 24 14:32:31 2012 + * + * @brief Test 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 "aka_common.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "material.hh" +#include "io_helper.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + initialize(argc, argv); + + // debug::setDebugLevel(dblDump); + + const UInt spatial_dimension = 2; + const UInt max_steps = 500; + + const ElementType type = _triangle_3; + + Mesh mesh(spatial_dimension); + MeshIOMSH mesh_io; + mesh_io.read("mesh.msh", mesh); + + SolidMechanicsModelCohesive model(mesh); + + /// model initialization + model.initExtrinsic("material.dat"); + Real time_step = model.getStableTimeStep()*0.05; + model.setTimeStep(time_step); + // std::cout << "Time step: " << time_step << std::endl; + + model.assembleMassLumped(); + + + /* ------------------------------------------------------------------------ */ + /* Facet part */ + /* ------------------------------------------------------------------------ */ + + // std::cout << mesh << std::endl; + + const Mesh & mesh_facets = model.getMeshFacets(); + + const ElementType type_facet = mesh.getFacetElementType(type); + UInt nb_facet = mesh_facets.getNbElement(type_facet); + // const Vector & position = mesh.getNodes(); + // const Vector & connectivity = mesh_facets.getConnectivity(type_facet); + + Vector & sigma_lim = model.getSigmaLimit(); + + Real * bary_facet = new Real[spatial_dimension]; + for (UInt f = 0; f < nb_facet; ++f) { + mesh_facets.getBarycenter(f, type_facet, bary_facet); + if (bary_facet[0] == -0.25) sigma_lim(f) = 100; + else sigma_lim(f) = 1e4; + } + delete[] bary_facet; + + // std::cout << mesh << std::endl; + + /* ------------------------------------------------------------------------ */ + /* End of facet part */ + /* ------------------------------------------------------------------------ */ + + // Vector & boundary = model.getBoundary(); + // const Vector & residual = model.getResidual(); + + UInt nb_nodes = mesh.getNbNodes(); + UInt nb_element = mesh.getNbElement(type); + + // /// boundary conditions + // for (UInt dim = 0; dim < spatial_dimension; ++dim) { + // for (UInt n = 0; n < nb_nodes; ++n) { + // boundary(n, dim) = true; + // } + // } + + model.updateResidual(); + + iohelper::ElemType paraview_type = iohelper::TRIANGLE1; + + /// initialize the paraview output + iohelper::DumperParaview dumper; + dumper.SetMode(iohelper::TEXT); + dumper.SetPoints(mesh.getNodes().values, + spatial_dimension, mesh.getNbNodes(), "explicit"); + dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, + paraview_type, nb_element, iohelper::C_MODE); + dumper.AddNodeDataField(model.getDisplacement().values, + spatial_dimension, "displacements"); + dumper.AddNodeDataField(model.getVelocity().values, + spatial_dimension, "velocity"); + dumper.AddNodeDataField(model.getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddNodeDataField(model.getForce().values, + spatial_dimension, "applied_force"); + dumper.AddNodeDataField(model.getResidual().values, + spatial_dimension, "forces"); + dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + spatial_dimension*spatial_dimension, "strain"); + dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + spatial_dimension*spatial_dimension, "stress"); + dumper.SetEmbeddedValue("displacements", 1); + dumper.SetEmbeddedValue("applied_force", 1); + dumper.SetEmbeddedValue("forces", 1); + dumper.SetPrefix("paraview/"); + dumper.Init(); + dumper.Dump(); + + + /// initial conditions + const Vector & position = mesh.getNodes(); + Vector & velocity = model.getVelocity(); + Real loading_rate = 1; + for (UInt n = 0; n < nb_nodes; ++n) { + velocity(n, 0) = loading_rate * position(n, 0); + } + + // std::ofstream edis("edis.txt"); + // std::ofstream erev("erev.txt"); + + Vector & residual = model.getResidual(); + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + + model.checkCohesiveStress(); + + model.explicitPred(); + model.updateResidual(); + model.updateAcceleration(); + model.explicitCorr(); + + if(s % 1 == 0) { + dumper.SetPoints(mesh.getNodes().values, + spatial_dimension, mesh.getNbNodes(), "explicit"); + dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, + paraview_type, nb_element, iohelper::C_MODE); + dumper.AddNodeDataField(model.getDisplacement().values, + spatial_dimension, "displacements"); + dumper.AddNodeDataField(model.getVelocity().values, + spatial_dimension, "velocity"); + dumper.AddNodeDataField(model.getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddNodeDataField(model.getForce().values, + spatial_dimension, "applied_force"); + dumper.AddNodeDataField(model.getResidual().values, + spatial_dimension, "forces"); + dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + spatial_dimension*spatial_dimension, "strain"); + dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + spatial_dimension*spatial_dimension, "stress"); + dumper.SetEmbeddedValue("displacements", 1); + dumper.SetEmbeddedValue("applied_force", 1); + dumper.SetEmbeddedValue("forces", 1); + dumper.SetPrefix("paraview/"); + dumper.Init(); + dumper.Dump(); + + std::cout << "passing step " << s << "/" << max_steps << std::endl; + } + + + // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); + + // edis << s << " " + // << Ed << std::endl; + + // erev << s << " " + // << Er << std::endl; + + } + + // edis.close(); + // erev.close(); + + Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + + Real Edt = 2; + + std::cout << Ed << " " << Edt << std::endl; + + // if (Ed < Edt - 0.0001 || Ed > Edt + 0.0001) { + // std::cout << "The dissipated energy is incorrect" << std::endl; + // return EXIT_FAILURE; + // } + + finalize(); + + std::cout << "OK: test_cohesive_extrinsic was passed!" << std::endl; + return EXIT_SUCCESS; +} diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/CMakeLists.txt new file mode 100644 index 000000000..da6dc01be --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/CMakeLists.txt @@ -0,0 +1,39 @@ +#=============================================================================== +# @file CMakeLists.txt +# @author Marco Vocialta +# @date Thu Apr 26 15:14:00 2012 +# +# @section LICENSE +# +# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) +# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) +# +# Akantu is free software: you can redistribute it and/or modify it under the +# terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR +# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Akantu. If not, see . +# +# @section DESCRIPTION +# +#=============================================================================== + +register_test(test_cohesive_intrinsic test_cohesive_intrinsic.cc) + +add_mesh(test_cohesive_intrinsic_mesh mesh.geo 2 2) +add_dependencies(test_cohesive_intrinsic test_cohesive_intrinsic_mesh) + +#=============================================================================== +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/material.dat + ${CMAKE_CURRENT_BINARY_DIR}/material.dat + COPYONLY + ) +#file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat new file mode 100644 index 000000000..54a145690 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/material.dat @@ -0,0 +1,15 @@ +material elastic [ + name = steel + rho = 1 # density + E = 1 # young's modulus + nu = 0.1 # poisson's ratio +] + +material cohesive_bilinear [ + name = cohesive + sigma_c = 1 + beta = 0 + G_cI = 1 + G_cII = 1 + delta_0 = 0.1 +] diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo new file mode 100644 index 000000000..4526cbf34 --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo @@ -0,0 +1,22 @@ +h = 1.0; + +Point(1) = {0, 0, 0, h}; +Point(2) = {0, 1, 0, h}; +Point(3) = {-1, 1, 0, h}; +Point(4) = {-1, 0, 0, h}; +Point(6) = {1, 0, 0, h}; +Point(7) = {1, 1, 0, h}; +Point(8) = {1, -1, 0, h}; +Point(9) = {0, -1, 0, h}; +Point(10) = {-1, -1, 0, h}; +Line(1) = {10, 9}; +Line(2) = {9, 8}; +Line(3) = {8, 6}; +Line(4) = {6, 7}; +Line(5) = {7, 2}; +Line(6) = {2, 3}; +Line(7) = {3, 4}; +Line(8) = {4, 10}; +Line Loop(9) = {3, 4, 5, 6, 7, 8, 1, 2}; +Plane Surface(10) = {9}; +Physical Surface(11) = {10}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc new file mode 100644 index 000000000..39906b52e --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc @@ -0,0 +1,240 @@ +/** + * @file test_cohesive.cc + * @author Marco Vocialta + * @date Fri Feb 24 14:32:31 2012 + * + * @brief Test 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 "aka_common.hh" +#include "mesh.hh" +#include "mesh_io.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "solid_mechanics_model_cohesive.hh" +#include "material.hh" +#include "io_helper.hh" +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +static void updateDisplacement(SolidMechanicsModelCohesive &, + Vector &, + ElementType, + Real); + +int main(int argc, char *argv[]) { + initialize(argc, argv); + + // debug::setDebugLevel(dblDump); + + const UInt spatial_dimension = 2; + const UInt max_steps = 350; + + const ElementType type = _triangle_6; + + Mesh mesh(spatial_dimension); + MeshIOMSH mesh_io; + mesh_io.read("mesh.msh", mesh); + + SolidMechanicsModelCohesive model(mesh); + + /// model initialization + model.initFull("material.dat"); + Real time_step = model.getStableTimeStep()*0.8; + model.setTimeStep(time_step); + // std::cout << "Time step: " << time_step << std::endl; + + model.assembleMassLumped(); + + + /* ------------------------------------------------------------------------ */ + /* Facet part */ + /* ------------------------------------------------------------------------ */ + + // std::cout << mesh << std::endl; + + const Mesh & mesh_facets = model.getMeshFacets(); + + const ElementType type_facet = mesh.getFacetElementType(type); + UInt nb_facet = mesh_facets.getNbElement(type_facet); + // const Vector & position = mesh.getNodes(); + // const Vector & connectivity = mesh_facets.getConnectivity(type_facet); + + Vector facet_insertion; + Real * bary_facet = new Real[spatial_dimension]; + for (UInt f = 0; f < nb_facet; ++f) { + mesh_facets.getBarycenter(f, type_facet, bary_facet); + if (bary_facet[0] == -0.25) facet_insertion.push_back(f); + } + delete[] bary_facet; + + model.insertCohesiveElements(facet_insertion); + + // mesh_io.write("mesh_cohesive.msh", mesh); + + // std::cout << mesh << std::endl; + + /* ------------------------------------------------------------------------ */ + /* End of facet part */ + /* ------------------------------------------------------------------------ */ + + Vector & boundary = model.getBoundary(); + // const Vector & residual = model.getResidual(); + + UInt nb_nodes = mesh.getNbNodes(); + UInt nb_element = mesh.getNbElement(type); + + /// boundary conditions + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + for (UInt n = 0; n < nb_nodes; ++n) { + boundary(n, dim) = true; + } + } + + model.updateResidual(); + + // iohelper::ElemType paraview_type = iohelper::TRIANGLE2; + + // /// initialize the paraview output + // iohelper::DumperParaview dumper; + // dumper.SetPoints(model.getFEM().getMesh().getNodes().values, + // spatial_dimension, nb_nodes, "explicit"); + // dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values, + // paraview_type, nb_element, iohelper::C_MODE); + // dumper.AddNodeDataField(model.getDisplacement().values, + // spatial_dimension, "displacements"); + // dumper.AddNodeDataField(model.getVelocity().values, + // spatial_dimension, "velocity"); + // dumper.AddNodeDataField(model.getAcceleration().values, + // spatial_dimension, "acceleration"); + // dumper.AddNodeDataField(model.getForce().values, + // spatial_dimension, "applied_force"); + // dumper.AddNodeDataField(model.getResidual().values, + // spatial_dimension, "forces"); + // dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + // spatial_dimension*spatial_dimension, "strain"); + // dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + // spatial_dimension*spatial_dimension, "stress"); + // dumper.SetEmbeddedValue("displacements", 1); + // dumper.SetEmbeddedValue("applied_force", 1); + // dumper.SetEmbeddedValue("forces", 1); + // dumper.SetPrefix("paraview/"); + // dumper.Init(); + // dumper.Dump(); + + + /// update displacement + Vector elements; + Real * bary = new Real[spatial_dimension]; + for (UInt el = 0; el < nb_element; ++el) { + mesh.getBarycenter(el, type, bary); + if (bary[0] > -0.25) elements.push_back(el); + } + delete[] bary; + + Real increment = 0.01; + + updateDisplacement(model, elements, type, increment); + + // std::ofstream edis("edis.txt"); + // std::ofstream erev("erev.txt"); + + /// Main loop + for (UInt s = 1; s <= max_steps; ++s) { + + model.explicitPred(); + model.updateResidual(); + model.updateAcceleration(); + model.explicitCorr(); + + if(s % 1 == 0) { + // dumper.Dump(); + std::cout << "passing step " << s << "/" << max_steps << std::endl; + } + + /// update displacement + updateDisplacement(model, elements, type, increment); + + // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); + + // edis << s << " " + // << Ed << std::endl; + + // erev << s << " " + // << Er << std::endl; + + } + + // edis.close(); + // erev.close(); + + Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); + + Real Edt = 2 * sqrt(2); + + if (Ed < Edt - 0.0001 || Ed > Edt + 0.0001) { + std::cout << "The dissipated energy is incorrect" << std::endl; + return EXIT_FAILURE; + } + + finalize(); + + std::cout << "OK: test_cohesive_intrinsic was passed!" << std::endl; + return EXIT_SUCCESS; +} + + +static void updateDisplacement(SolidMechanicsModelCohesive & model, + Vector & elements, + ElementType type, + Real increment) { + + Mesh & mesh = model.getFEM().getMesh(); + UInt nb_element = elements.getSize(); + UInt nb_nodes = mesh.getNbNodes(); + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + + const Vector & connectivity = mesh.getConnectivity(type); + Vector & displacement = model.getDisplacement(); + Vector update(nb_nodes); + + for (UInt el = 0; el < nb_element; ++el) { + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt node = connectivity(elements(el), n); + if (!update(node)) { + displacement(node, 0) += increment; + // displacement(node, 1) -= increment; + update(node) = true; + } + } + } +}