Diffusion akantu (features/node_to_node) Browse

node_to_node

material_igfem_iterative_stiffness_reduction.cc
View Options

			/**
			* @file material_igfem_iterative_stiffness_reduction.cc
			* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
			* @date Thu Mar 10 08:37:43 2016
			*
			* @brief Implementation of igfem material iterative stiffness reduction
			*
			*
			* 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 <http://www.gnu.org/licenses/>.
			*
			*/

			/* -------------------------------------------------------------------------- */
			#include "material_igfem_iterative_stiffness_reduction.hh"
			#include <math.h>

			namespace akantu {
			template <Int spatial_dimension>

			/* -------------------------------------------------------------------------- */
			MaterialIGFEMIterativeStiffnessReduction<spatial_dimension>::
			MaterialIGFEMIterativeStiffnessReduction(SolidMechanicsModel & model,
			const ID & id)
			: Material(model, id), MaterialIGFEMSawToothDamage<spatial_dimension>(model,
			id),
			eps_u("ultimate_strain", this), reduction_step("damage_step", this),
			D("tangent", *this), Gf(0.), crack_band_width(0.), max_reductions(0),
			reduction_constant(0.) {
			AKANTU_DEBUG_IN();

			this->eps_u.initialize(1);
			this->D.initialize(1);
			this->reduction_step.initialize(1);

			this->internals_to_transfer.push_back("ultimate_strain");
			this->internals_to_transfer.push_back("tangent");

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <Int dim>
			void MaterialIGFEMIterativeStiffnessReduction<dim>::initMaterial() {
			AKANTU_DEBUG_IN();

			MaterialIGFEMSawToothDamage<dim>::initMaterial();

			/// get the parameters for the sub-material that can be damaged
			ID mat_name = this->sub_material_names[1];
			const Material & mat = this->model->getMaterial(mat_name);
			this->crack_band_width = mat.getParam<Real>("crack_band_width");
			this->max_reductions = mat.getParam<UInt>("max_reductions");
			this->reduction_constant = mat.getParam<Real>("reduction_constant");
			this->Gf = mat.getParam<Real>("Gf");

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <Int spatial_dimension>
			void MaterialIGFEMIterativeStiffnessReduction<spatial_dimension>::
			computeNormalizedEquivalentStress(const Array<Real> & grad_u,
			ElementType el_type,
			GhostType ghost_type) {
			AKANTU_DEBUG_IN();
			/// storage for the current stress
			Matrix<Real> sigma(spatial_dimension, spatial_dimension);
			/// Vector to store eigenvalues of current stress tensor
			Vector<Real> eigenvalues(spatial_dimension);

			/// iterators on the needed internal fields
			Array<Real>::const_scalar_iterator Sc_it =
			this->Sc(el_type, ghost_type).begin();
			Array<Real>::scalar_iterator dam_it =
			this->damage(el_type, ghost_type).begin();
			Array<Real>::scalar_iterator equivalent_stress_it =
			this->equivalent_stress(el_type, ghost_type).begin();
			Array<Real>::const_matrix_iterator grad_u_it =
			grad_u.begin(spatial_dimension, spatial_dimension);
			Array<Real>::const_matrix_iterator grad_u_end =
			grad_u.end(spatial_dimension, spatial_dimension);
			Real * mu_ptr = this->mu(el_type, ghost_type).data();
			Real * lambda_ptr = this->lambda(el_type, ghost_type).data();

			/// loop over all the quadrature points and compute the equivalent stress
			for (; grad_u_it != grad_u_end; ++grad_u_it) {
			/// compute the stress
			sigma.zero();
			MaterialIGFEMElastic<spatial_dimension>::computeStressOnQuad(
			grad_u_it, sigma, lambda_ptr, *mu_ptr);
			MaterialIGFEMSawToothDamage<
			spatial_dimension>::computeDamageAndStressOnQuad(sigma, *dam_it);
			/// compute eigenvalues
			sigma.eig(eigenvalues);

			/// find max eigenvalue and normalize by tensile strength
			*equivalent_stress_it =
			*(std::max_element(eigenvalues.data(),
			eigenvalues.data() + spatial_dimension)) /
			(*Sc_it);
			++Sc_it;
			++equivalent_stress_it;
			++dam_it;
			++lambda_ptr;
			++mu_ptr;
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <Int spatial_dimension>
			UInt MaterialIGFEMIterativeStiffnessReduction<
			spatial_dimension>::updateDamage() {
			UInt nb_damaged_elements = 0;

			if (this->norm_max_equivalent_stress >= 1.) {

			AKANTU_DEBUG_IN();

			/// update the damage only on non-ghosts elements! Doesn't make sense to
			/// update on ghost.
			GhostType ghost_type = _not_ghost;
			;

			Mesh::type_iterator it = this->model->getFEEngine().getMesh().firstType(
			spatial_dimension, ghost_type, _ek_igfem);
			Mesh::type_iterator last_type =
			this->model->getFEEngine().getMesh().lastType(spatial_dimension,
			ghost_type, _ek_igfem);

			/// get the Young's modulus of the damageable sub-material
			ID mat_name = this->sub_material_names[1];

			Real E = this->model->getMaterial(mat_name).template getParam<Real>("E");

			/// loop over all the elements
			for (; it != last_type; ++it) {
			ElementType el_type = *it;

			/// get iterators on the needed internal fields
			const Array<Idx> & sub_mat = this->sub_material(el_type, ghost_type);
			Array<Idx>::const_scalar_iterator sub_mat_it = sub_mat.begin();
			Array<Real>::const_scalar_iterator equivalent_stress_it =
			this->equivalent_stress(el_type, ghost_type).begin();
			Array<Real>::const_scalar_iterator equivalent_stress_end =
			this->equivalent_stress(el_type, ghost_type).end();
			Array<Real>::scalar_iterator dam_it =
			this->damage(el_type, ghost_type).begin();
			Array<Idx>::scalar_iterator reduction_it =
			this->reduction_step(el_type, ghost_type).begin();
			Array<Real>::const_scalar_iterator eps_u_it =
			this->eps_u(el_type, ghost_type).begin();
			Array<Real>::scalar_iterator Sc_it =
			this->Sc(el_type, ghost_type).begin();
			Array<Real>::const_scalar_iterator D_it =
			this->D(el_type, ghost_type).begin();

			/// loop over all the elements of the given type
			UInt nb_element = this->element_filter(el_type, ghost_type).getSize();
			UInt nb_quads = this->fem->getNbIntegrationPoints(el_type, ghost_type);
			bool damage_element = false;
			for (Int e = 0; e < nb_element; ++e) {
			damage_element = false;
			/// check if damage occurs in the element
			for (Int q = 0; q < nb_quads;
			++q, ++reduction_it, ++sub_mat_it, ++equivalent_stress_it) {
			if (equivalent_stress_it >= (1 - this->dam_tolerance)
			this->norm_max_equivalent_stress &&
			*sub_mat_it != 0) {
			/// check if this element can still be damaged
			if (*reduction_it == this->max_reductions)
			continue;
			damage_element = true;
			}
			}

			if (damage_element) {
			/// damage the element
			nb_damaged_elements += 1;
			sub_mat_it -= nb_quads;
			reduction_it -= nb_quads;
			for (Int q = 0; q < nb_quads; ++q) {
			if (*sub_mat_it) {
			/// increment the counter of stiffness reduction steps
			*reduction_it += 1;
			if (*reduction_it == this->max_reductions)
			*dam_it = this->max_damage;
			else {
			/// update the damage on this quad
			*dam_it = 1. - (1. / std::pow(this->reduction_constant,
			*reduction_it));
			/// update the stiffness on this quad
			Sc_it = (eps_u_it) * (1. - (dam_it)) E * (*D_it) /
			((1. - (dam_it)) E + (*D_it));
			}
			}
			++sub_mat_it;
			++dam_it;
			++reduction_it;
			++eps_u_it;
			++Sc_it;
			++D_it;
			}
			} else {
			dam_it += nb_quads;
			eps_u_it += nb_quads;
			Sc_it += nb_quads;
			D_it += nb_quads;
			}
			}
			}
			}
			StaticCommunicator & comm =
			akantu::StaticCommunicator::getStaticCommunicator();
			comm.allReduce(&nb_damaged_elements, 1, _so_sum);
			AKANTU_DEBUG_OUT();
			return nb_damaged_elements;
			}

			/* -------------------------------------------------------------------------- */
			template <Int spatial_dimension>
			void MaterialIGFEMIterativeStiffnessReduction<
			spatial_dimension>::onElementsAdded(__attribute__((unused))
			const Array<Element> & element_list,
			__attribute__((unused))
			const NewElementsEvent & event) {

			MaterialIGFEMSawToothDamage<spatial_dimension>::onElementsAdded(element_list,
			event);
			/// set the correct damage iteration step (is UInt->cannot be interpolated)
			Real val = 0.;
			for (ghost_type_t::iterator g = ghost_type_t::begin();
			g != ghost_type_t::end(); ++g) {
			GhostType ghost_type = *g;
			/// loop over all types in the material
			typedef ElementTypeMapArray<Idx>::type_iterator iterator;
			iterator it = this->element_filter.firstType(spatial_dimension, ghost_type,
			_ek_igfem);
			iterator last_type =
			this->element_filter.lastType(spatial_dimension, ghost_type, _ek_igfem);
			/// loop over all types in the filter
			for (; it != last_type; ++it) {
			const ElementType el_type = *it;
			Array<Real>::scalar_iterator dam_it =
			this->damage(el_type, ghost_type).begin();
			Array<Idx>::scalar_iterator reduction_it =
			this->reduction_step(el_type, ghost_type).begin();
			UInt nb_element = this->element_filter(el_type, ghost_type).getSize();
			UInt nb_quads = this->fem->getNbIntegrationPoints(el_type);
			UInt * sub_mat_ptr = this->sub_material(el_type, ghost_type).data();
			for (Int q = 0; q < nb_element * nb_quads;
			++q, ++sub_mat_ptr, ++dam_it, ++reduction_it) {
			if (*sub_mat_ptr) {
			if (Math::are_float_equal(*dam_it, this->max_damage))
			*reduction_it = this->max_reductions;
			else {
			for (Int i = 0; i < this->max_reductions; ++i) {
			val = 1 - (1. / std::pow(this->reduction_constant, i));
			if (Math::are_float_equal(val, *dam_it))
			*reduction_it = i;
			}
			}
			}
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */

			INSTANTIATE_MATERIAL(MaterialIGFEMIterativeStiffnessReduction);

			} // namespace akantu

rAKA/extra_packages/igfem/src/material_igfem/material_igfem_iterative_stiffness_reduction.cce29c5d5d51cbfeatures/node_to_node

material_igfem_iterative_stiffness_reduction.ccView Options

rAKA/extra_packages/igfem/src/material_igfem/material_igfem_iterative_stiffness_reduction.cc
e29c5d5d51cbfeatures/node_to_node

material_igfem_iterative_stiffness_reduction.cc
View Options