diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc index 1b6fd6c7f..36432dff2 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc @@ -1,363 +1,367 @@ /** * @file material_cohesive_exponential.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Mon Jul 09 2012 * @date last modification: Tue Aug 04 2015 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_exponential.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveExponential::MaterialCohesiveExponential(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id) { AKANTU_DEBUG_IN(); - this->registerParam("beta" , beta , Real(0.) , _pat_parsable, "Beta parameter"); + this->registerParam("beta" , beta , Real(0.), + _pat_parsable, + "Beta parameter"); + this->registerParam("exponential_penalty", exp_penalty, true, - _pat_parsable, "Is contact penalty following the exponential law?" ); + _pat_parsable, + "Is contact penalty following the exponential law?"); this->registerParam("contact_tangent", contact_tangent, Real(1.0), - _pat_parsable, "Ratio of contact tangent over the initial exponential tangent" ); + _pat_parsable, + "Ratio of contact tangent over the initial exponential tangent"); // this->initInternalArray(delta_max, 1, _ek_cohesive); use_previous_delta_max=true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); - if ((exp_penalty)&&(contact_tangent!=1)){ + if ((exp_penalty) && (contact_tangent != 1)){ contact_tangent = 1; - AKANTU_DEBUG_WARNING("The parsed paramter is forced to 1.0 when the contact penalty follows the exponential law"); + AKANTU_DEBUG_WARNING("The parsed paramter is forced to 1.0 when the +contact penalty follows the exponential law"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator traction_end = tractions(el_type, ghost_type).end(spatial_dimension); Array::iterator delta_max_it = delta_max(el_type, ghost_type).begin(); Array::iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); /// compute scalars - Real beta2 = beta*beta; + Real beta2 = beta * beta; /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++delta_max_it, ++delta_max_prev_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \beta^2 \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); - if ((normal_opening_norm<0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { + if ((normal_opening_norm < 0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { - Vector op_n(*opening_it); - op_n *= *normal_it; + Vector op_n(*normal_it); + op_n *= normal_opening_norm; Vector delta_s(*opening_it); delta_s -= op_n; - delta = tangential_opening_norm*beta; + delta = tangential_opening_norm * beta; computeCoupledTraction(*traction_it, *normal_it, delta, delta_s, *delta_max_it, *delta_max_prev_it); computeCompressiveTraction(*traction_it, *normal_it, normal_opening_norm, *opening_it); } else computeCoupledTraction(*traction_it, *normal_it, delta, *opening_it, *delta_max_it, *delta_max_prev_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTraction(Vector & tract, const Vector & normal, Real delta, const Vector & opening, Real & delta_max_new, Real delta_max) { AKANTU_DEBUG_IN(); /// full damage case if (std::abs(delta) < Math::getTolerance()) { /// set traction to zero tract.clear(); } else { /// element not fully damaged /** * Compute traction loading @f$ \mathbf{T} = * e \sigma_c \frac{\delta}{\delta_c} e^{-\delta/ \delta_c}@f$ */ /** * Compute traction unloading @f$ \mathbf{T} = * \frac{t_{max}}{\delta_{max}} \delta @f$ */ - Real beta2 = beta*beta; - Vector op_n_n(opening); - op_n_n*=normal; - op_n_n*=(1-beta2); - op_n_n*=normal; - tract = beta2*opening; + Real beta2 = beta * beta; + Real normal_open_norm = opening.dot(normal); + Vector op_n_n(spatial_dimension); + op_n_n = normal; + op_n_n *= (1 - beta2); + op_n_n *= normal_open_norm; + tract = beta2 * opening; tract += op_n_n; delta_max_new = std::max(delta_max, delta); - tract *= exp(1)*sigma_c*exp(- delta_max_new / delta_c)/delta_c; + tract *= exp(1) * sigma_c * exp(-delta_max_new / delta_c) / delta_c; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressiveTraction(Vector & tract, const Vector & normal, Real delta_n, __attribute__((unused)) const Vector & opening) { Vector temp_tract(normal); if(exp_penalty) { - temp_tract *= delta_n*exp(1)*sigma_c*exp(- delta_n / delta_c)/delta_c; + temp_tract *= delta_n * exp(1) * sigma_c * exp(-delta_n / delta_c) / delta_c; } else { - Real initial_tg = contact_tangent*exp(1)*sigma_c*delta_n/delta_c; + Real initial_tg = contact_tangent * exp(1) * sigma_c * delta_n / delta_c; temp_tract *= initial_tg; } tract += temp_tract; } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array::matrix_iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); - Array::matrix_iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); - Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); - Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); - Array::iteratordelta_max_it = delta_max.previous(el_type, ghost_type).begin(); - Real beta2 = beta*beta; + Array::matrix_iterator tangent_it = + tangent_matrix.begin(spatial_dimension, spatial_dimension); + + Array::matrix_iterator tangent_end = + tangent_matrix.end(spatial_dimension, spatial_dimension); + + Array::const_vector_iterator normal_it = + normal.begin(spatial_dimension); + + Array::vector_iterator opening_it = + opening(el_type, ghost_type).begin(spatial_dimension); + + Array::iteratordelta_max_it = + delta_max.previous(el_type, ghost_type).begin(); + + Real beta2 = beta * beta; /** * compute tangent matrix @f$ \frac{\partial \mathbf{t}} * {\partial \delta} = \hat{\mathbf{t}} \otimes * \frac{\partial (t/\delta)}{\partial \delta} * \frac{\hat{\mathbf{t}}}{\delta}+ \frac{t}{\delta} [ \beta^2 \mathbf{I} + * (1-\beta^2) (\mathbf{n} \otimes \mathbf{n})] @f$ **/ /** * In which @f$ * \frac{\partial(t/ \delta)}{\partial \delta} = * \left\{\begin{array} {l l} * -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad if * \delta \geq \delta_{max} \\ * 0 & \quad if \delta < \delta_{max}, \delta_n > 0 * \end{array}\right. @f$ **/ - for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++traction_it, ++ delta_max_it) { + for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, + ++opening_it, ++ delta_max_it) { + Real normal_opening_norm = opening_it->dot(*normal_it); Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; - Real tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); - if ((normal_opening_norm<0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { + if ((normal_opening_norm < 0) && (std::abs(normal_opening_norm) > Math::getTolerance())) { - Vector op_n(*opening_it); - op_n *= *normal_it; + Vector op_n(*normal_it); + op_n *= normal_opening_norm; Vector delta_s(*opening_it); delta_s -= op_n; - delta = tangential_opening_norm*beta; + delta = tangential_opening_norm * beta; - computeCoupledTangent(*tangent_it, *traction_it, *normal_it, delta, - delta_s, *delta_max_it); + computeCoupledTangent(*tangent_it, *normal_it, delta, delta_s, *delta_max_it); + computeCompressivePenalty(*tangent_it, *normal_it, normal_opening_norm); } else computeCoupledTangent(*tangent_it, *traction_it, *normal_it, delta, *opening_it, *delta_max_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCoupledTangent(Matrix & tangent, - Vector tract, const Vector & normal, Real delta, const Vector & opening, Real delta_max_new) { AKANTU_DEBUG_IN(); - Real beta2 = beta*beta; - Matrix I(spatial_dimension, spatial_dimension); - I.eye(beta2); - Matrix nn(spatial_dimension, spatial_dimension); - nn.outerProduct(normal, normal); - - nn *= (1-beta2); - I += nn; + Real beta2 = beta * beta; + Matrix J(spatial_dimension,spatial_dimension); + J.eye(beta2); if(std::abs(delta) < Math::getTolerance()){ - I *= exp(1)*sigma_c/delta_c; - tangent += I; - } else { - - Real traction_norm = tract.norm(); - - Vector t_hat(opening); - Vector beta2delta(opening); - - beta2delta *= beta2; - t_hat *= normal; - t_hat *= (beta2-1); - t_hat *= normal; - t_hat += beta2delta; - - Vector deriv_t_hat(t_hat); - - deriv_t_hat /= delta; - - Real derivatives = 0; - - if ((delta > delta_max_new) || (std::abs(delta - delta_max_new) < 1e-12)) { - derivatives = -exp(1- delta/delta_c)*sigma_c/(delta_c*delta_c); + delta = Math::getTolerance(); } + + Real opening_normal; + opening_normal = opening.dot(normal); - deriv_t_hat *= derivatives; - - Matrix t_var_t(spatial_dimension, spatial_dimension); - t_var_t.outerProduct(t_hat, deriv_t_hat); - - I *= traction_norm / delta; + Vectordelta_e(normal); + delta_e *= opening_normal; + delta_e *= (1 - beta2); + delta_e += (beta2 * opening); + + Real exponent = exp(1 - delta / delta_c) * sigma_c / delta_c; + + Matrix first_term(spatial_dimension, spatial_dimension); + first_term.outerProduct(normal, normal); + first_term *= (1 - beta2); + first_term += J; + + Matrixsecond_term(spatial_dimension, spatial_dimension); + second_term.outerProduct(delta_e,delta_e); + second_term /= delta; + second_term /= delta_c; - tangent += t_var_t; - tangent += I; + Matrixdiff(first_term); + diff -= second_term; - } + tangent = diff; + tangent *= exponent; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeCompressivePenalty(Matrix & tangent, const Vector & normal, Real delta_n) { if (!exp_penalty) delta_n=0; Matrix n_outer_n(spatial_dimension,spatial_dimension); n_outer_n.outerProduct(normal,normal); - Real normal_tg = contact_tangent*exp(1)*sigma_c*exp(-delta_n/delta_c)*(1-delta_n/delta_c)/delta_c; + Real normal_tg = contact_tangent * exp(1) * sigma_c * exp(-delta_n / delta_c) * (1 - delta_n / delta_c) / delta_c; n_outer_n *= normal_tg; tangent += n_outer_n; } INSTANTIATE_MATERIAL(MaterialCohesiveExponential); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh index 942c138ed..11de18644 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh @@ -1,128 +1,128 @@ /** * @file material_cohesive_exponential.hh * * @author Fabian Barras * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 06 2015 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ #define __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material Exponential damage * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: 0) * - delta_c : critical opening (default: 0) */ template class MaterialCohesiveExponential : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveExponential(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// Initialization void initMaterial(); /// constitutive law void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type = _not_ghost); private: void computeCoupledTraction(Vector & tract, const Vector & normal, Real delta, const Vector & opening, Real & delta_max_new, Real delta_max); void computeCompressiveTraction(Vector & tract, const Vector & normal, Real delta_n, const Vector & opening); - void computeCoupledTangent(Matrix & tangent, Vector tract, + void computeCoupledTangent(Matrix & tangent, const Vector & normal, Real delta, const Vector & opening, Real delta_max_new); void computeCompressivePenalty(Matrix & tangent, const Vector & normal, Real delta_n); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// beta parameter Real beta; /// contact penalty = initial slope ? bool exp_penalty; /// Ratio of contact tangent over the initial exponential tangent Real contact_tangent; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ // #include "material_cohesive_exponential_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc index a56f19514..bc07e92f2 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,625 +1,625 @@ /** * @file material_cohesive.cc * * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jan 12 2016 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" #include "aka_random_generator.hh" #include "shape_cohesive.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive(&(model.getFEEngineClass("CohesiveFEEngine"))), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), opening_old("opening (old)", *this), tractions("tractions", *this), tractions_old("tractions (old)", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta max", *this), use_previous_delta_max(false), use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normal(0, spatial_dimension, "normal") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->model->getMesh().initElementTypeMapArray(this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, spatial_dimension - 1); this->reversible_energy.initialize(1 ); this->total_energy .initialize(1 ); this->tractions_old .initialize(spatial_dimension); this->tractions .initialize(spatial_dimension); this->opening_old .initialize(spatial_dimension); this->contact_tractions.initialize(spatial_dimension); this->contact_opening .initialize(spatial_dimension); this->opening .initialize(spatial_dimension); this->delta_max .initialize(1 ); this->damage .initialize(1 ); if (this->model->getIsExtrinsic()) this->sigma_c.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) this->delta_max.initializeHistory(); if (this->use_previous_opening) this->opening.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif Array & residual = const_cast &>(model->getResidual()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & traction = tractions(*it, ghost_type); Array & contact_traction = contact_tractions(*it, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem_cohesive->getNbIntegrationPoints(*it, ghost_type); /// compute @f$t_i N_a@f$ Array * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); Array::iterator > traction_it = traction.begin(spatial_dimension, 1); Array::iterator > contact_traction_it = contact_traction.begin(spatial_dimension, 1); Array::const_iterator > shapes_filtered_begin = shapes.begin(1, size_of_shapes); Array::iterator > traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ Array * int_t_N = new Array(nb_element, spatial_dimension*size_of_shapes, "int_t_N"); fem_cohesive->integrate(*traction_cpy, *int_t_N, spatial_dimension * size_of_shapes, *it, ghost_type, elem_filter); delete traction_cpy; int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent()); Real * int_t_N_val = int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element * spatial_dimension; } /// assemble model->getFEEngineBoundary().assembleArray(*int_t_N, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, 1); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { UInt nb_quadrature_points = fem_cohesive->getNbIntegrationPoints(*it, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if(!nb_element) continue; UInt size_of_shapes = shapes.getNbComponent(); Array * shapes_filtered = new Array(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_val = shapes.storage(); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } /** * compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c c c} * 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ * 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ * 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ * 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ * 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\ * 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 * \end{array} \right]@f$ **/ // UInt size_of_A = spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element; // Real * A = new Real[size_of_A]; // memset(A, 0, size_of_A*sizeof(Real)); Matrix A(spatial_dimension*size_of_shapes, spatial_dimension*nb_nodes_per_element); for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension*size_of_shapes) = -1; } - /// compute traction. This call is not necessary for the linear - /// cohesive law that, currently, is the only one used for the - /// extrinsic approach. - if (!model->getIsExtrinsic()){ - computeTraction(ghost_type); - } + // compute traction. This call is not necessary for the linear + // cohesive law that, currently, is the only one used for the + // extrinsic approach. + // if (!model->getIsExtrinsic()){ + // computeTraction(ghost_type); + //} /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} @f$ Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * nb_quadrature_points, spatial_dimension, "normal"); normal.resize(nb_quadrature_points); computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ //computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); tangent_stiffness_matrix->clear(); computeTangentTraction(*it, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension*nb_nodes_per_element*spatial_dimension*nb_nodes_per_element; Array * at_nt_d_n_a = new Array (nb_element*nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator > shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N (spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A (spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for(; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.clear(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension ; ++i) for (UInt n = 0; n < size_of_shapes; ++n) N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} {\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; Array * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive->integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, *it, ghost_type, elem_filter); delete at_nt_d_n_a; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, *it, ghost_type, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = nb_element*fem_cohesive->getNbIntegrationPoints(*it, ghost_type); normal.resize(nb_quadrature_points); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (type == _cohesive_1d_2) fem_cohesive->computeNormalsOnIntegrationPoints(position, normal, type, ghost_type); else { #define COMPUTE_NORMAL(type) \ fem_cohesive->getShapeFunctions(). \ computeNormalsOnIntegrationPoints(position, \ normal, \ ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_OPENING(type) \ fem_cohesive->getShapeFunctions(). \ interpolateOnIntegrationPoints(displacement, \ opening, \ spatial_dimension, \ ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); Real * memory_space = new Real[2*spatial_dimension]; Vector b(memory_space, spatial_dimension); Vector h(memory_space + spatial_dimension, spatial_dimension); for(; it != last_type; ++it) { Array::iterator erev = reversible_energy(*it, _not_ghost).begin(); Array::iterator etot = total_energy(*it, _not_ghost).begin(); Array::vector_iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_end = tractions(*it, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } delete [] memory_space; /// update old values it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); GhostType ghost_type = _not_ghost; for(; it != last_type; ++it) { tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type)); opening_old(*it, ghost_type).copy(opening(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Array dissipated_energy(total_energy(*it, _not_ghost)); dissipated_energy -= reversible_energy(*it, _not_ghost); edis += fem_cohesive->integrate(dissipated_energy, *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Array & el_filter = element_filter(*it, _not_ghost); UInt nb_quad_per_el = fem_cohesive->getNbIntegrationPoints(*it, _not_ghost); UInt nb_quad_points = el_filter.getSize() * nb_quad_per_el; Array contact_energy(nb_quad_points); Array::vector_iterator contact_traction_it = contact_tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(*it, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt el = 0; el < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++el) { contact_energy(el) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive->integrate(contact_energy, *it, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); else if (type == "cohesive contact") return getContactEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ __END_AKANTU__