solid_mechanics_model_RVE.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Nov 28, 14:23

solid_mechanics_model_RVE.cc
View Options

	/**
	* @file solid_mechanics_model_RVE.cc
	* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
	* @date Wed Jan 13 15:32:35 2016
	*
	* @brief Implementation of SolidMechanicsModelRVE
	*
	* @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 <http://www.gnu.org/licenses/>.
	*
	*/

	/* -------------------------------------------------------------------------- */
	#include "solid_mechanics_model_RVE.hh"
	#include "material_elastic.hh"
	#ifdef AKANTU_USE_MUMPS
	#include "solver_mumps.hh"
	#endif
	/* -------------------------------------------------------------------------- */
	__BEGIN_AKANTU__

	/* -------------------------------------------------------------------------- */
	SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh, bool use_RVE_mat_selector, UInt dim, const ID & id, const MemoryID & memory_id) :
	SolidMechanicsModel(mesh, dim, id, memory_id),
	volume(0.),
	use_RVE_mat_selector(use_RVE_mat_selector),
	static_communicator_dummy(StaticCommunicator::getStaticCommunicatorDummy()) {

	/// create node groups for PBCs
	mesh.createGroupsFromMeshData<std::string>("physical_names");
	/// find the four corner nodes of the RVE
	findCornerNodes();

	/// remove the corner nodes from the surface node groups:
	/// This most be done because corner nodes a not periodic
	mesh.getElementGroup("top").removeNode(corner_nodes(2));
	mesh.getElementGroup("top").removeNode(corner_nodes(3));
	mesh.getElementGroup("left").removeNode(corner_nodes(3));
	mesh.getElementGroup("left").removeNode(corner_nodes(0));
	mesh.getElementGroup("bottom").removeNode(corner_nodes(1));
	mesh.getElementGroup("bottom").removeNode(corner_nodes(0));
	mesh.getElementGroup("right").removeNode(corner_nodes(2));
	mesh.getElementGroup("right").removeNode(corner_nodes(1));

	const ElementGroup & bottom = mesh.getElementGroup("bottom");
	bottom_nodes.insert( bottom.node_begin(), bottom.node_end() );

	const ElementGroup & left = mesh.getElementGroup("left");
	left_nodes.insert( left.node_begin(), left.node_end() );

	/// enforce periodicity on the displacement fluctuations
	SurfacePair surface_pair_1 = std::make_pair("top","bottom");
	SurfacePair surface_pair_2 = std::make_pair("right","left");
	SurfacePairList surface_pairs_list;
	surface_pairs_list.push_back(surface_pair_1);
	surface_pairs_list.push_back(surface_pair_2);
	this->setPBC(surface_pairs_list);
	}

	/* -------------------------------------------------------------------------- */
	SolidMechanicsModelRVE::~SolidMechanicsModelRVE() {
	delete static_communicator_dummy;
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::initFull(const ModelOptions & options) {
	SolidMechanicsModel::initFull(options);

	this->initMaterials();

	/// compute the volume of the RVE
	FEEngine * fem = this->fems["SolidMechanicsFEEngine"];
	GhostType gt = _not_ghost;
	Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
	Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
	for(; it != end; ++it) {
	const ElementType element_type = *it;
	Array<Real> Volume(this->mesh.getNbElement(element_type) * this->fems["SolidMechanicsFEEngine"]->getNbIntegrationPoints(element_type), 1, 1.);
	this->volume = fem->integrate(Volume, element_type);
	}

	std::cout << "The volume of the RVE is " << this->volume << std::endl;

	/// dumping
	std::stringstream base_name;
	base_name << "RVE_" << this->memory_id - 1;
	this->setBaseName (base_name.str());
	this->addDumpFieldVector("displacement");
	this->addDumpField ("stress" );
	this->addDumpField ("grad_u" );
	this->addDumpField ("blocked_dofs" );
	this->addDumpField ("material_index" );

	this->dump();
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::applyBoundaryConditions(const Matrix<Real> & displacement_gradient) {

	/// get the position of the nodes
	const Array<Real> & pos = mesh.getNodes();
	/// storage for the coordinates of a given node and the displacement that will be applied
	Vector<Real> x(spatial_dimension);
	Vector<Real> appl_disp(spatial_dimension);

	/// fix top right node
	UInt node = this->corner_nodes(2);
	x(0) = pos(node,0); x(1) = pos(node,1);
	appl_disp.mul<false>(displacement_gradient,x);
	(this->blocked_dofs)(node,0) = true; (this->displacement)(node,0) = appl_disp(0);
	(this->blocked_dofs)(node,1) = true; (this->displacement)(node,1) = appl_disp(1);
	// (this->blocked_dofs)(node,0) = true; (this->displacement)(node,0) = 0.;
	// (this->blocked_dofs)(node,1) = true; (this->displacement)(node,1) = 0.;

	/// apply Hx at all the other corner nodes; H: displ. gradient
	node = this->corner_nodes(0);
	x(0) = pos(node,0); x(1) = pos(node,1);
	appl_disp.mul<false>(displacement_gradient,x);
	(this->blocked_dofs)(node,0) = true; (this->displacement)(node,0) = appl_disp(0);
	(this->blocked_dofs)(node,1) = true; (this->displacement)(node,1) = appl_disp(1);

	node = this->corner_nodes(1);
	x(0) = pos(node,0); x(1) = pos(node,1);
	appl_disp.mul<false>(displacement_gradient,x);
	(this->blocked_dofs)(node,0) = true; (this->displacement)(node,0) = appl_disp(0);
	(this->blocked_dofs)(node,1) = true; (this->displacement)(node,1) = appl_disp(1);

	node = this->corner_nodes(3);
	x(0) = pos(node,0); x(1) = pos(node,1);
	appl_disp.mul<false>(displacement_gradient,x);
	(this->blocked_dofs)(node,0) = true; (this->displacement)(node,0) = appl_disp(0);
	(this->blocked_dofs)(node,1) = true; (this->displacement)(node,1) = appl_disp(1);
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::findCornerNodes() {
	mesh.computeBoundingBox();

	// find corner nodes
	const Array<Real> & position = mesh.getNodes();
	const Vector<Real> & lower_bounds = mesh.getLowerBounds();
	const Vector<Real> & upper_bounds = mesh.getUpperBounds();

	AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!");
	corner_nodes.resize(4);
	corner_nodes.set(UInt(-1));

	for (UInt n = 0; n < mesh.getNbNodes(); ++n) {
	// node 1
	if (Math::are_float_equal(position(n, 0), lower_bounds(0)) &&
	Math::are_float_equal(position(n, 1), lower_bounds(1))) {
	corner_nodes(0) = n;
	}
	// node 2
	else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) &&
	Math::are_float_equal(position(n, 1), lower_bounds(1))) {
	corner_nodes(1) = n;
	}
	// node 3
	else if (Math::are_float_equal(position(n, 0), upper_bounds(0)) &&
	Math::are_float_equal(position(n, 1), upper_bounds(1))) {
	corner_nodes(2) = n;
	}
	// node 4
	else if (Math::are_float_equal(position(n, 0), lower_bounds(0)) &&
	Math::are_float_equal(position(n, 1), upper_bounds(1))) {
	corner_nodes(3) = n;
	}
	}

	for (UInt i = 0; i < corner_nodes.getSize(); ++i) {
	if (corner_nodes(i) == UInt(-1))
	AKANTU_DEBUG_ERROR("The corner node " << i + 1 << " wasn't found");
	}
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::advanceASR(const Matrix<Real> & prestrain) {

	/// apply the new eigenstrain
	GhostType gt = _not_ghost;

	Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
	Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
	for(; it != end; ++it) {
	const ElementType element_type = *it;
	Array<Real> & prestrain_vect = const_cast<Array<Real> &>(this->getMaterial("gel").getInternal<Real>("eigen_grad_u")(element_type, gt));
	Array<Real>::iterator< Matrix<Real> > prestrain_it = prestrain_vect.begin(spatial_dimension, spatial_dimension);
	Array<Real>::iterator< Matrix<Real> > prestrain_end = prestrain_vect.end(spatial_dimension, spatial_dimension);

	for (; prestrain_it != prestrain_end; ++prestrain_it)
	(*prestrain_it) = prestrain;
	}

	/// solve the system for the given boundary conditions
	Real error = 0;
	bool converged = this->solveStep<_scm_newton_raphson_tangent, _scc_increment>(1e-10, error, 2, false, *static_communicator_dummy);
	AKANTU_DEBUG_ASSERT(converged, "Did not converge");
	}


	/* -------------------------------------------------------------------------- */
	Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index, const ID & field_type) {

	FEEngine * fem = this->fems["SolidMechanicsFEEngine"];
	Real average = 0;
	GhostType gt = _not_ghost;

	Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
	Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
	for(; it != end; ++it) {
	const ElementType element_type = *it;
	if (field_type == "stress") {
	for(UInt m = 0; m < this->materials.size(); ++m) {
	const Array<Real> & stress_vec = this->materials[m]->getStress(element_type);
	const Array<UInt> & elem_filter = this->materials[m]->getElementFilter(element_type);
	Array<Real> int_stress_vec(elem_filter.getSize(), spatial_dimension*spatial_dimension, "int_of_stress");

	fem->integrate(stress_vec, int_stress_vec, spatial_dimension*spatial_dimension, element_type, _not_ghost, elem_filter);

	for (UInt k = 0; k < elem_filter.getSize(); ++k)
	average += int_stress_vec(k, row_index * spatial_dimension + col_index); //3 is the value for the yy (in 3D, the value is 4)
	}
	}

	else if (field_type == "strain") {
	for(UInt m = 0; m < this->materials.size(); ++m) {
	const Array<Real> & gradu_vec = this->materials[m]->getGradU(element_type);
	const Array<UInt> & elem_filter = this->materials[m]->getElementFilter(element_type);
	Array<Real> int_gradu_vec(elem_filter.getSize(), spatial_dimension*spatial_dimension, "int_of_gradu");

	fem->integrate(gradu_vec, int_gradu_vec, spatial_dimension*spatial_dimension, element_type, _not_ghost, elem_filter);

	for (UInt k = 0; k < elem_filter.getSize(); ++k)
	/// averaging is done only for normal components, so stress and strain are equal
	average += 0.5 * (int_gradu_vec(k, row_index * spatial_dimension + col_index) + int_gradu_vec(k, col_index * spatial_dimension + row_index));
	}
	}

	else if (field_type == "eigen_grad_u") {
	for(UInt m = 0; m < this->materials.size(); ++m) {
	const Array<Real> & eigen_gradu_vec = this->materials[m]->getInternal<Real>("eigen_grad_u")(element_type);
	const Array<UInt> & elem_filter = this->materials[m]->getElementFilter(element_type);
	Array<Real> int_eigen_gradu_vec(elem_filter.getSize(), spatial_dimension*spatial_dimension, "int_of_gradu");

	fem->integrate(eigen_gradu_vec, int_eigen_gradu_vec, spatial_dimension*spatial_dimension, element_type, _not_ghost, elem_filter);

	for (UInt k = 0; k < elem_filter.getSize(); ++k)
	/// averaging is done only for normal components, so stress and strain are equal
	average += int_eigen_gradu_vec(k, row_index * spatial_dimension + col_index);
	}
	}

	else
	AKANTU_DEBUG_ERROR("Averaging not implemented for this field!!!");
	}
	return average/this->volume;
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::homogenizeStiffness(Matrix<Real> & C_macro) {

	const UInt dim = 2;
	AKANTU_DEBUG_ASSERT(this->spatial_dimension == dim, "Is only implemented for 2D!!!");

	/// apply three independent loading states to determine C
	/// 1. eps_el = (1;0;0) 2. eps_el = (0,1,0) 3. eps_el = (0,0,0.5)

	/// storage for results of 3 different loading states
	UInt voigt_size = VoigtHelper<dim>::size;
	Matrix<Real> stresses(voigt_size, voigt_size, 0.);
	Matrix<Real> strains(voigt_size, voigt_size, 0.);
	Matrix<Real> H(dim, dim, 0.);

	/// virtual test 1:
	H(0,0) = 1;
	this->performVirtualTesting(H, stresses, strains, 0);

	/// virtual test 2:
	H.clear();
	H(1,1) = 1.;
	this->performVirtualTesting(H, stresses, strains, 1);

	/// virtual test 3:
	H.clear();
	H(0,1) = 1.;
	this->performVirtualTesting(H, stresses, strains, 2);

	/// compute effective stiffness
	Matrix<Real> eps_inverse(voigt_size, voigt_size);
	eps_inverse.inverse(strains);
	C_macro.mul<false, false>(stresses, eps_inverse);
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::performVirtualTesting(const Matrix<Real> & H, Matrix<Real> & eff_stresses, Matrix<Real> & eff_strains, const UInt test_no) {

	this->applyBoundaryConditions(H);

	/// solve system
	this->assembleStiffnessMatrix();
	Real error = 0;
	bool converged= this->solveStep<_scm_newton_raphson_tangent_not_computed, _scc_increment>(1e-12, error, 2, false, *static_communicator_dummy);
	AKANTU_DEBUG_ASSERT(converged, "Did not converge");

	/// get average stress and strain
	eff_stresses(0, test_no) = this->averageTensorField(0,0, "stress");
	eff_strains(0, test_no) = this->averageTensorField(0,0, "strain");
	eff_stresses(1, test_no) = this->averageTensorField(1,1, "stress");
	eff_strains(1, test_no) = this->averageTensorField(1,1, "strain");
	eff_stresses(2, test_no) = this->averageTensorField(1,0, "stress");
	eff_strains(2, test_no) = 2. * this->averageTensorField(1,0, "strain");

	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::homogenizeEigenGradU(Matrix<Real> & eigen_gradu_macro) {
	eigen_gradu_macro(0,0) = this->averageTensorField(0,0, "eigen_grad_u");
	eigen_gradu_macro(1,1) = this->averageTensorField(1,1, "eigen_grad_u");
	eigen_gradu_macro(0,1) = this->averageTensorField(0,1, "eigen_grad_u");
	eigen_gradu_macro(1,0) = this->averageTensorField(1,0, "eigen_grad_u");
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::initMaterials() {
	AKANTU_DEBUG_IN();

	// make sure the material are instantiated
	if(!are_materials_instantiated) instantiateMaterials();

	if (use_RVE_mat_selector) {
	const Vector<Real> & lowerBounds = mesh.getLowerBounds();
	const Vector<Real> & upperBounds = mesh.getUpperBounds();
	Real bottom = lowerBounds(1);
	Real top = upperBounds(1);
	Real box_size = std::abs(top-bottom);
	Real eps = box_size * 1e-6;


	if(is_default_material_selector) delete material_selector;
	material_selector = new GelMaterialSelector(*this, box_size, "gel", 1, eps);
	is_default_material_selector = false;
	}

	SolidMechanicsModel::initMaterials();

	AKANTU_DEBUG_OUT();
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::initSolver(__attribute__((unused)) SolverOptions & options) {
	#if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
	AKANTU_DEBUG_ERROR("You should at least activate one solver.");
	#else
	UInt nb_global_nodes = mesh.getNbGlobalNodes();

	delete jacobian_matrix;
	std::stringstream sstr; sstr << id << ":jacobian_matrix";

	#ifdef AKANTU_USE_PETSC
	jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
	#else
	jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _unsymmetric, sstr.str(), memory_id, 1);
	#endif //AKANTU_USE PETSC
	jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);

	if (!isExplicit()) {
	delete stiffness_matrix;
	std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
	#ifdef AKANTU_USE_PETSC
	stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id, 1);
	stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
	#else
	stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
	#endif //AKANTU_USE_PETSC
	}

	delete solver;
	std::stringstream sstr_solv; sstr_solv << id << ":solver";
	#ifdef AKANTU_USE_PETSC
	solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str(), memory_id);
	#elif defined(AKANTU_USE_MUMPS)
	solver = new SolverMumps(*jacobian_matrix, sstr_solv.str(), memory_id,
	*static_communicator_dummy);
	dof_synchronizer->initScatterGatherCommunicationScheme();
	#else
	AKANTU_DEBUG_ERROR("You should at least activate one solver.");
	#endif //AKANTU_USE_MUMPS

	SolverMumpsOptions opt(SolverMumpsOptions::_not_parallel);

	if(solver)
	solver->initialize(opt);
	#endif //AKANTU_HAS_SOLVER
	}

	/* -------------------------------------------------------------------------- */
	void SolidMechanicsModelRVE::initArrays() {
	AKANTU_DEBUG_IN();

	UInt nb_nodes = mesh.getNbNodes();
	std::stringstream sstr_disp; sstr_disp << id << ":displacement";
	// std::stringstream sstr_mass; sstr_mass << id << ":mass";
	std::stringstream sstr_velo; sstr_velo << id << ":velocity";
	std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
	std::stringstream sstr_forc; sstr_forc << id << ":force";
	std::stringstream sstr_resi; sstr_resi << id << ":residual";
	std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs";

	displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
	// mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, spatial_dimension, 0));
	velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
	acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
	force = &(alloc<Real>(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
	residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
	blocked_dofs = &(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));

	std::stringstream sstr_curp; sstr_curp << id << ":current_position";
	current_position = &(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));

	for(UInt g = _not_ghost; g <= _ghost; ++g) {
	GhostType gt = (GhostType) g;
	Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
	Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
	for(; it != end; ++it) {
	UInt nb_element = mesh.getNbElement(*it, gt);
	material_index.alloc(nb_element, 1, *it, gt);
	material_local_numbering.alloc(nb_element, 1, *it, gt);
	}
	}

	dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension,
	*static_communicator_dummy);
	dof_synchronizer->initLocalDOFEquationNumbers();
	dof_synchronizer->initGlobalDOFEquationNumbers();

	AKANTU_DEBUG_OUT();
	}


	__END_AKANTU__

solid_mechanics_model_RVE.ccNo OneTemporaryActions

File Metadata

solid_mechanics_model_RVE.ccView Options

Event Timeline

solid_mechanics_model_RVE.cc
No OneTemporary
Actions

solid_mechanics_model_RVE.cc
View Options