Diffusion akantu (bugfixes/46-bug-in-strain-energy-split-and-missing-constraint-in-phasefield) Browse

46-bug-in…

fragment_manager.cc
View Options

			/**
			* @file fragment_manager.cc
			*
			* @author Marco Vocialta <marco.vocialta@epfl.ch>
			*
			* @date creation: Thu Jan 23 2014
			* @date last modification: Tue Aug 19 2014
			*
			* @brief Group manager to handle fragments
			*
			* @section LICENSE
			*
			* Copyright (©) 2014 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 "fragment_manager.hh"
			#include "material_cohesive.hh"
			#include <numeric>
			#include <algorithm>
			#include <functional>

			__BEGIN_AKANTU__

			/* -------------------------------------------------------------------------- */
			FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model,
			bool dump_data,
			const ID & id,
			const MemoryID & memory_id) :
			GroupManager(model.getMesh(), id, memory_id),
			model(model),
			mass_center(0, model.getSpatialDimension(), "mass_center"),
			mass(0, model.getSpatialDimension(), "mass"),
			velocity(0, model.getSpatialDimension(), "velocity"),
			inertia_moments(0, model.getSpatialDimension(), "inertia_moments"),
			principal_directions(0, model.getSpatialDimension() * model.getSpatialDimension(),
			"principal_directions"),
			quad_coordinates("quad_coordinates", id),
			mass_density("mass_density", id),
			nb_elements_per_fragment(0, 1, "nb_elements_per_fragment"),
			dump_data(dump_data) {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = mesh.getSpatialDimension();

			/// compute quadrature points' coordinates
			mesh.initElementTypeMapArray(quad_coordinates,
			spatial_dimension,
			spatial_dimension,
			_not_ghost);

			model.getFEEngine().interpolateOnQuadraturePoints(model.getMesh().getNodes(),
			quad_coordinates);

			/// store mass density per quadrature point
			mesh.initElementTypeMapArray(mass_density,
			1,
			spatial_dimension,
			_not_ghost);

			storeMassDensityPerQuadraturePoint();

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			class CohesiveElementFilter : public GroupManager::ClusteringFilter {
			public:
			CohesiveElementFilter(const SolidMechanicsModelCohesive & model,
			const Real max_damage = 1.) :
			model(model), is_unbroken(max_damage) {}

			bool operator() (const Element & el) const {
			if (el.kind == _ek_regular)
			return true;

			const Array<UInt> & mat_indexes = model.getMaterialByElement(el.type,
			el.ghost_type);
			const Array<UInt> & mat_loc_num = model.getMaterialLocalNumbering(el.type,
			el.ghost_type);

			const MaterialCohesive & mat
			= static_cast<const MaterialCohesive &>
			(model.getMaterial(mat_indexes(el.element)));

			UInt el_index = mat_loc_num(el.element);
			UInt nb_quad_per_element
			= model.getFEEngine("CohesiveFEEngine").getNbQuadraturePoints(el.type, el.ghost_type);

			const Array<Real> & damage_array = mat.getDamage(el.type, el.ghost_type);

			AKANTU_DEBUG_ASSERT(nb_quad_per_element * el_index < damage_array.getSize(),
			"This quadrature point is out of range");

			const Real * element_damage
			= damage_array.storage() + nb_quad_per_element * el_index;

			UInt unbroken_quads = std::count_if(element_damage,
			element_damage + nb_quad_per_element,
			is_unbroken);

			if (unbroken_quads > 0)
			return true;
			return false;
			}

			private:

			struct IsUnbrokenFunctor {
			IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {}
			bool operator() (const Real & x) {return x < max_damage;}
			const Real max_damage;
			};

			const SolidMechanicsModelCohesive & model;
			const IsUnbrokenFunctor is_unbroken;
			};

			/* -------------------------------------------------------------------------- */
			void FragmentManager::buildFragments(Real damage_limit) {
			AKANTU_DEBUG_IN();

			#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
			DistributedSynchronizer * cohesive_synchronizer
			= const_cast<DistributedSynchronizer *>(model.getCohesiveSynchronizer());

			if (cohesive_synchronizer) {
			cohesive_synchronizer->computeBufferSize(model, _gst_smmc_damage);
			cohesive_synchronizer->asynchronousSynchronize(model, _gst_smmc_damage);
			cohesive_synchronizer->waitEndSynchronize(model, _gst_smmc_damage);
			}
			#endif

			DistributedSynchronizer & synchronizer
			= const_cast<DistributedSynchronizer &>(model.getSynchronizer());

			Mesh & mesh_facets = const_cast<Mesh &>(mesh.getMeshFacets());

			UInt spatial_dimension = model.getSpatialDimension();
			std::string fragment_prefix("fragment");

			/// generate fragments
			global_nb_fragment = createClusters(spatial_dimension,
			fragment_prefix,
			CohesiveElementFilter(model,
			damage_limit),
			&synchronizer,
			&mesh_facets);

			nb_fragment = getNbElementGroups(spatial_dimension);
			fragment_index.resize(nb_fragment);

			UInt * fragment_index_it = fragment_index.storage();

			/// loop over fragments
			for(const_element_group_iterator it(element_group_begin());
			it != element_group_end(); ++it, ++fragment_index_it) {

			/// get fragment index
			std::string fragment_index_string
			= it->first.substr(fragment_prefix.size() + 1);
			std::stringstream sstr(fragment_index_string.c_str());
			sstr >> *fragment_index_it;

			AKANTU_DEBUG_ASSERT(!sstr.fail(), "fragment_index is not an integer");
			}

			/// compute fragments' mass
			computeMass();

			if (dump_data) {
			createDumpDataArray(fragment_index, "fragments", true);
			createDumpDataArray(mass, "fragments mass");
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::computeMass() {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = model.getSpatialDimension();

			/// create a unit field per quadrature point, since to compute mass
			/// it's neccessary to integrate only density
			ElementTypeMapArray<Real> unit_field("unit_field", id);
			mesh.initElementTypeMapArray(unit_field,
			spatial_dimension,
			spatial_dimension,
			_not_ghost);

			ElementTypeMapArray<Real>::type_iterator it = unit_field.firstType(spatial_dimension,
			_not_ghost,
			_ek_regular);
			ElementTypeMapArray<Real>::type_iterator end = unit_field.lastType(spatial_dimension,
			_not_ghost,
			_ek_regular);

			for (; it != end; ++it) {
			ElementType type = *it;
			Array<Real> & field_array = unit_field(type);
			UInt nb_element = mesh.getNbElement(type);
			UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type);

			field_array.resize(nb_element * nb_quad_per_element);
			field_array.set(1.);
			}

			integrateFieldOnFragments(unit_field, mass);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::computeCenterOfMass() {
			AKANTU_DEBUG_IN();

			/// integrate position multiplied by density
			integrateFieldOnFragments(quad_coordinates, mass_center);

			/// divide it by the fragments' mass
			Real * mass_storage = mass.storage();
			Real * mass_center_storage = mass_center.storage();

			UInt total_components = mass_center.getSize() * mass_center.getNbComponent();

			for (UInt i = 0; i < total_components; ++i)
			mass_center_storage[i] /= mass_storage[i];

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::computeVelocity() {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = model.getSpatialDimension();

			/// compute velocity per quadrature point
			ElementTypeMapArray<Real> velocity_field("velocity_field", id);

			mesh.initElementTypeMapArray(velocity_field,
			spatial_dimension,
			spatial_dimension,
			_not_ghost);

			model.getFEEngine().interpolateOnQuadraturePoints(model.getVelocity(),
			velocity_field);

			/// integrate on fragments
			integrateFieldOnFragments(velocity_field, velocity);

			/// divide it by the fragments' mass
			Real * mass_storage = mass.storage();
			Real * velocity_storage = velocity.storage();

			UInt total_components = velocity.getSize() * velocity.getNbComponent();

			for (UInt i = 0; i < total_components; ++i)
			velocity_storage[i] /= mass_storage[i];

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			/**
			* Given the distance @f$ \mathbf{r} @f$ between a quadrature point
			* and its center of mass, the moment of inertia is computed as \f[
			* I_\mathrm{CM} = \mathrm{tr}(\mathbf{r}\mathbf{r}^\mathrm{T})
			* \mathbf{I} - \mathbf{r}\mathbf{r}^\mathrm{T} \f] for more
			* information check Wikipedia
			* (http://en.wikipedia.org/wiki/Moment_of_inertia#Identities_for_a_skew-symmetric_matrix)
			*
			*/

			void FragmentManager::computeInertiaMoments() {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = model.getSpatialDimension();

			computeCenterOfMass();

			/// compute local coordinates products with respect to the center of match
			ElementTypeMapArray<Real> moments_coords("moments_coords", id);

			mesh.initElementTypeMapArray(moments_coords,
			spatial_dimension * spatial_dimension,
			spatial_dimension,
			_not_ghost);

			/// resize the by element type
			ElementTypeMapArray<Real>::type_iterator it = moments_coords.firstType(spatial_dimension,
			_not_ghost,
			_ek_regular);
			ElementTypeMapArray<Real>::type_iterator end = moments_coords.lastType(spatial_dimension,
			_not_ghost,
			_ek_regular);

			for (; it != end; ++it) {
			ElementType type = *it;
			Array<Real> & field_array = moments_coords(type);
			UInt nb_element = mesh.getNbElement(type);
			UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type);

			field_array.resize(nb_element * nb_quad_per_element);
			}


			/// compute coordinates
			Array<Real>::const_vector_iterator mass_center_it
			= mass_center.begin(spatial_dimension);

			/// loop over fragments
			for(const_element_group_iterator it(element_group_begin());
			it != element_group_end(); ++it, ++mass_center_it) {

			const ElementTypeMapArray<UInt> & el_list = it->second->getElements();

			ElementTypeMapArray<UInt>::type_iterator type_it = el_list.firstType(spatial_dimension,
			_not_ghost,
			_ek_regular);
			ElementTypeMapArray<UInt>::type_iterator type_end = el_list.lastType(spatial_dimension,
			_not_ghost,
			_ek_regular);

			/// loop over elements of the fragment
			for (; type_it != type_end; ++type_it) {
			ElementType type = *type_it;
			UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type);

			Array<Real> & moments_coords_array = moments_coords(type);
			const Array<Real> & quad_coordinates_array = quad_coordinates(type);
			const Array<UInt> & el_list_array = el_list(type);

			Array<Real>::matrix_iterator moments_begin
			= moments_coords_array.begin(spatial_dimension, spatial_dimension);
			Array<Real>::const_vector_iterator quad_coordinates_begin
			= quad_coordinates_array.begin(spatial_dimension);

			Vector<Real> relative_coords(spatial_dimension);

			for (UInt el = 0; el < el_list_array.getSize(); ++el) {
			UInt global_el = el_list_array(el);

			/// loop over quadrature points
			for (UInt q = 0; q < nb_quad_per_element; ++q) {
			UInt global_q = global_el * nb_quad_per_element + q;
			Matrix<Real> moments_matrix = moments_begin[global_q];
			const Vector<Real> & quad_coord_vector = quad_coordinates_begin[global_q];

			/// to understand this read the documentation written just
			/// before this function
			relative_coords = quad_coord_vector;
			relative_coords -= *mass_center_it;

			moments_matrix.outerProduct(relative_coords, relative_coords);
			Real trace = moments_matrix.trace();
			moments_matrix *= -1.;
			moments_matrix += Matrix<Real>::eye(spatial_dimension, trace);
			}
			}
			}
			}

			/// integrate moments
			Array<Real> integrated_moments(global_nb_fragment,
			spatial_dimension * spatial_dimension);

			integrateFieldOnFragments(moments_coords, integrated_moments);

			/// compute and store principal moments
			inertia_moments.resize(global_nb_fragment);
			principal_directions.resize(global_nb_fragment);

			Array<Real>::matrix_iterator integrated_moments_it
			= integrated_moments.begin(spatial_dimension, spatial_dimension);
			Array<Real>::vector_iterator inertia_moments_it
			= inertia_moments.begin(spatial_dimension);
			Array<Real>::matrix_iterator principal_directions_it
			= principal_directions.begin(spatial_dimension, spatial_dimension);

			for (UInt frag = 0; frag < global_nb_fragment; ++frag, ++integrated_moments_it,
			++inertia_moments_it, ++principal_directions_it) {
			integrated_moments_it->eig(inertia_moments_it, principal_directions_it);
			}

			if (dump_data)
			createDumpDataArray(inertia_moments, "moments of inertia");

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::computeAllData() {
			AKANTU_DEBUG_IN();

			buildFragments();
			computeVelocity();
			computeInertiaMoments();

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::storeMassDensityPerQuadraturePoint() {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = model.getSpatialDimension();

			Mesh::type_iterator it = mesh.firstType(spatial_dimension);
			Mesh::type_iterator end = mesh.lastType(spatial_dimension);

			for (; it != end; ++it) {
			ElementType type = *it;

			Array<Real> & mass_density_array = mass_density(type);

			UInt nb_element = mesh.getNbElement(type);
			UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type);
			mass_density_array.resize(nb_element * nb_quad_per_element);

			const Array<UInt> & mat_indexes = model.getMaterialByElement(type);

			Real * mass_density_it = mass_density_array.storage();

			/// store mass_density for each element and quadrature point
			for (UInt el = 0; el < nb_element; ++el) {
			Material & mat = model.getMaterial(mat_indexes(el));

			for (UInt q = 0; q < nb_quad_per_element; ++q, ++mass_density_it)
			*mass_density_it = mat.getRho();
			}
			}

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::integrateFieldOnFragments(ElementTypeMapArray<Real> & field,
			Array<Real> & output) {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = model.getSpatialDimension();
			UInt nb_component = output.getNbComponent();

			/// integration part
			output.resize(global_nb_fragment);
			output.clear();

			UInt * fragment_index_it = fragment_index.storage();
			Array<Real>::vector_iterator output_begin = output.begin(nb_component);

			/// loop over fragments
			for(const_element_group_iterator it(element_group_begin());
			it != element_group_end(); ++it, ++fragment_index_it) {

			const ElementTypeMapArray<UInt> & el_list = it->second->getElements();

			ElementTypeMapArray<UInt>::type_iterator type_it = el_list.firstType(spatial_dimension,
			_not_ghost,
			_ek_regular);
			ElementTypeMapArray<UInt>::type_iterator type_end = el_list.lastType(spatial_dimension,
			_not_ghost,
			_ek_regular);

			/// loop over elements of the fragment
			for (; type_it != type_end; ++type_it) {
			ElementType type = *type_it;
			UInt nb_quad_per_element = model.getFEEngine().getNbQuadraturePoints(type);

			const Array<Real> & density_array = mass_density(type);
			Array<Real> & field_array = field(type);
			const Array<UInt> & elements = el_list(type);
			UInt nb_element = elements.getSize();

			/// generate array to be integrated by filtering fragment's elements
			Array<Real> integration_array(elements.getSize() * nb_quad_per_element,
			nb_component);

			Array<Real>::matrix_iterator int_array_it
			= integration_array.begin_reinterpret(nb_quad_per_element,
			nb_component, nb_element);
			Array<Real>::matrix_iterator int_array_end
			= integration_array.end_reinterpret(nb_quad_per_element,
			nb_component, nb_element);
			Array<Real>::matrix_iterator field_array_begin
			= field_array.begin_reinterpret(nb_quad_per_element,
			nb_component,
			field_array.getSize() / nb_quad_per_element);
			Array<Real>::const_vector_iterator density_array_begin
			= density_array.begin_reinterpret(nb_quad_per_element,
			density_array.getSize() / nb_quad_per_element);

			for (UInt el = 0; int_array_it != int_array_end; ++int_array_it, ++el) {
			UInt global_el = elements(el);
			*int_array_it = field_array_begin[global_el];

			/// multiply field by density
			const Vector<Real> & density_vector = density_array_begin[global_el];

			for (UInt i = 0; i < nb_quad_per_element; ++i) {
			for (UInt j = 0; j < nb_component; ++j) {
			(int_array_it)(i, j) = density_vector(i);
			}
			}
			}


			/// integrate the field over the fragment
			Array<Real> integrated_array(elements.getSize(), nb_component);
			model.getFEEngine().integrate(integration_array,
			integrated_array,
			nb_component,
			type,
			_not_ghost,
			elements);

			/// sum over all elements and store the result
			Vector<Real> output_tmp(output_begin[*fragment_index_it]);
			output_tmp += std::accumulate(integrated_array.begin(nb_component),
			integrated_array.end(nb_component),
			Vector<Real>(nb_component));
			}
			}

			/// sum output over all processors
			StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
			comm.allReduce(output.storage(), global_nb_fragment * nb_component, _so_sum);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void FragmentManager::computeNbElementsPerFragment() {
			AKANTU_DEBUG_IN();

			UInt spatial_dimension = model.getSpatialDimension();
			nb_elements_per_fragment.resize(global_nb_fragment);
			nb_elements_per_fragment.clear();

			UInt * fragment_index_it = fragment_index.storage();

			/// loop over fragments
			for(const_element_group_iterator it(element_group_begin());
			it != element_group_end(); ++it, ++fragment_index_it) {

			const ElementTypeMapArray<UInt> & el_list = it->second->getElements();

			ElementTypeMapArray<UInt>::type_iterator type_it = el_list.firstType(spatial_dimension,
			_not_ghost,
			_ek_regular);
			ElementTypeMapArray<UInt>::type_iterator type_end = el_list.lastType(spatial_dimension,
			_not_ghost,
			_ek_regular);

			/// loop over elements of the fragment
			for (; type_it != type_end; ++type_it) {
			ElementType type = *type_it;
			UInt nb_element = el_list(type).getSize();

			nb_elements_per_fragment(*fragment_index_it) += nb_element;
			}
			}

			/// sum values over all processors
			StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
			comm.allReduce(nb_elements_per_fragment.storage(), global_nb_fragment, _so_sum);

			if (dump_data)
			createDumpDataArray(nb_elements_per_fragment, "elements per fragment");

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			template <typename T>
			void FragmentManager::createDumpDataArray(Array<T> & data,
			std::string name,
			bool fragment_index_output) {
			AKANTU_DEBUG_IN();

			if (data.getSize() == 0) return;

			typedef typename Array<T>::vector_iterator data_iterator;
			Mesh & mesh_not_const = const_cast<Mesh &>(mesh);

			UInt spatial_dimension = mesh.getSpatialDimension();
			UInt nb_component = data.getNbComponent();
			UInt * fragment_index_it = fragment_index.storage();

			data_iterator data_begin = data.begin(nb_component);

			/// loop over fragments
			for(const_element_group_iterator it(element_group_begin());
			it != element_group_end(); ++it, ++fragment_index_it) {

			const ElementGroup & fragment = *(it->second);

			/// loop over cluster types
			ElementGroup::type_iterator type_it = fragment.firstType(spatial_dimension);
			ElementGroup::type_iterator type_end = fragment.lastType(spatial_dimension);

			for(; type_it != type_end; ++type_it) {
			ElementType type = *type_it;

			/// init mesh data
			Array<T> * mesh_data
			= mesh_not_const.getDataPointer<T>(name, type, _not_ghost, nb_component);

			data_iterator mesh_data_begin = mesh_data->begin(nb_component);

			ElementGroup::const_element_iterator el_it = fragment.element_begin(type);
			ElementGroup::const_element_iterator el_end = fragment.element_end(type);

			/// fill mesh data
			if (fragment_index_output) {
			for (; el_it != el_end; ++el_it) {
			Vector<T> md_tmp(mesh_data_begin[*el_it]);
			md_tmp(0) = *fragment_index_it;
			}
			} else {
			for (; el_it != el_end; ++el_it) {
			Vector<T> md_tmp(mesh_data_begin[*el_it]);
			md_tmp = data_begin[*fragment_index_it];
			}
			}
			}
			}

			AKANTU_DEBUG_OUT();
			}


			__END_AKANTU__

rAKA/src/model/solid_mechanics/fragment_manager.cc9054f18b641ebugfixes/46-bug-in…

fragment_manager.ccView Options

rAKA/src/model/solid_mechanics/fragment_manager.cc
9054f18b641ebugfixes/46-bug-in…

fragment_manager.cc
View Options