Diffusion akantu (bugfixes/65-Making-the-integrators-available-in-Python) Browse

65-Making-the…

solid_mechanics_model_RVE.cc
View Options

			/**
			* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
			* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
			*
			* This file is part of Akantu
			*
			* 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 "element_group.hh"
			#include "material_damage_iterative.hh"
			#include "node_group.hh"
			#include "non_linear_solver.hh"
			#include "non_local_manager.hh"
			#include "parser.hh"
			#include "sparse_matrix.hh"
			/* -------------------------------------------------------------------------- */
			#include <string>
			/* -------------------------------------------------------------------------- */

			namespace akantu {

			/* -------------------------------------------------------------------------- */
			SolidMechanicsModelRVE::SolidMechanicsModelRVE(Mesh & mesh,
			bool use_RVE_mat_selector,
			UInt nb_gel_pockets, Int dim,
			const ID & id)
			: SolidMechanicsModel(mesh, dim, id), volume(0.),
			use_RVE_mat_selector(use_RVE_mat_selector),
			nb_gel_pockets(nb_gel_pockets), nb_dumps(0) {
			AKANTU_DEBUG_IN();
			/// 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 auto & bottom = mesh.getElementGroup("bottom").getNodeGroup();
			bottom_nodes.insert(bottom.begin(), bottom.end());

			const auto & left = mesh.getElementGroup("left").getNodeGroup();
			left_nodes.insert(left.begin(), left.end());

			// /// enforce periodicity on the displacement fluctuations
			// auto surface_pair_1 = std::make_pair("top", "bottom");
			// auto 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);
			// TODO: To Nicolas correct the PBCs
			// this->setPBC(surface_pairs_list);

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default;

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) {
			AKANTU_DEBUG_IN();

			auto options_cp(options);
			options_cp.analysis_method = AnalysisMethod::_static;

			SolidMechanicsModel::initFullImpl(options_cp);

			// this->initMaterials();
			auto & fem = this->getFEEngine("SolidMechanicsFEEngine");

			/// compute the volume of the RVE
			GhostType gt = _not_ghost;
			for (auto element_type :
			this->mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) {
			Array<Real> Volume(this->mesh.getNbElement(element_type) *
			fem.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 << this->id;
			this->setBaseName(base_name.str());
			this->addDumpFieldVector("displacement");
			this->addDumpField("stress");
			this->addDumpField("grad_u");
			this->addDumpField("eigen_grad_u");
			this->addDumpField("blocked_dofs");
			this->addDumpField("material_index");
			this->addDumpField("damage");
			this->addDumpField("Sc");
			this->addDumpField("external_force");
			this->addDumpField("equivalent_stress");
			this->addDumpField("internal_force");
			this->addDumpField("delta_T");

			this->dump();
			this->nb_dumps += 1;
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::applyBoundaryConditions(
			const Matrix<Real> & displacement_gradient) {
			AKANTU_DEBUG_IN();
			/// 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);
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::applyHomogeneousTemperature(
			const Real & temperature) {

			for (UInt m = 0; m < this->getNbMaterials(); ++m) {
			Material & mat = this->getMaterial(m);

			const ElementTypeMapArray<Idx> & filter_map = mat.getElementFilter();

			// Loop over all element types
			for (auto && type : filter_map.elementTypes(spatial_dimension)) {
			const Array<Idx> & filter = filter_map(type);
			if (filter.size() == 0)
			continue;

			Array<Real> & delta_T = mat.getArray<Real>("delta_T", type);
			Array<Real>::scalar_iterator delta_T_it = delta_T.begin();
			Array<Real>::scalar_iterator delta_T_end = delta_T.end();

			for (; delta_T_it != delta_T_end; ++delta_T_it) {
			*delta_T_it = temperature;
			}
			}
			}
			}

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

			// find corner nodes
			const auto & position = mesh.getNodes();
			const auto & lower_bounds = mesh.getLowerBounds();
			const auto & upper_bounds = mesh.getUpperBounds();

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

			for (auto && data : enumerate(make_view(position, spatial_dimension))) {
			auto node = std::get<0>(data);
			const auto & X = std::get<1>(data);

			auto distance = X.distance(lower_bounds);
			// node 1
			if (Math::are_float_equal(distance, 0)) {
			corner_nodes(0) = node;
			}
			// node 2
			else if (Math::are_float_equal(X(_x), upper_bounds(_x)) &&
			Math::are_float_equal(X(_y), lower_bounds(_y))) {
			corner_nodes(1) = node;
			}
			// node 3
			else if (Math::are_float_equal(X(_x), upper_bounds(_x)) &&
			Math::are_float_equal(X(_y), upper_bounds(_y))) {
			corner_nodes(2) = node;
			}
			// node 4
			else if (Math::are_float_equal(X(_x), lower_bounds(_x)) &&
			Math::are_float_equal(X(_y), upper_bounds(_y))) {
			corner_nodes(3) = node;
			}
			}

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

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::advanceASR(const Matrix<Real> & prestrain) {
			AKANTU_DEBUG_IN();
			AKANTU_DEBUG_ASSERT(spatial_dimension == 2, "This is 2D only!");

			/// apply the new eigenstrain
			for (auto element_type :
			mesh.elementTypes(spatial_dimension, _not_ghost, _ek_not_defined)) {
			Array<Real> & prestrain_vect =
			const_cast<Array<Real> &>(this->getMaterial("gel").getInternal<Real>(
			"eigen_grad_u")(element_type));
			auto prestrain_it =
			prestrain_vect.begin(spatial_dimension, spatial_dimension);
			auto prestrain_end =
			prestrain_vect.end(spatial_dimension, spatial_dimension);

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

			/// advance the damage
			MaterialDamageIterative<2> & mat_paste =
			dynamic_cast<MaterialDamageIterative<2> &>(*this->materials[1]);
			MaterialDamageIterative<2> & mat_aggregate =
			dynamic_cast<MaterialDamageIterative<2> &>(*this->materials[0]);
			UInt nb_damaged_elements = 0;
			Real max_eq_stress_aggregate = 0;
			Real max_eq_stress_paste = 0;

			auto & solver = this->getNonLinearSolver();
			solver.set("max_iterations", 2);
			solver.set("threshold", 1e-6);
			solver.set("convergence_type", SolveConvergenceCriteria::_solution);

			do {
			this->solveStep();

			/// compute damage
			max_eq_stress_aggregate = mat_aggregate.getNormMaxEquivalentStress();
			max_eq_stress_paste = mat_paste.getNormMaxEquivalentStress();

			nb_damaged_elements = 0;
			if (max_eq_stress_aggregate > max_eq_stress_paste)
			nb_damaged_elements = mat_aggregate.updateDamage();
			else if (max_eq_stress_aggregate < max_eq_stress_paste)
			nb_damaged_elements = mat_paste.updateDamage();
			else
			nb_damaged_elements =
			(mat_paste.updateDamage() + mat_aggregate.updateDamage());

			std::cout << "the number of damaged elements is " << nb_damaged_elements
			<< std::endl;
			} while (nb_damaged_elements);

			if (this->nb_dumps % 10 == 0) {
			this->dump();
			}
			this->nb_dumps += 1;

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			Real SolidMechanicsModelRVE::averageTensorField(UInt row_index, UInt col_index,
			const ID & field_type) {
			AKANTU_DEBUG_IN();
			auto & fem = this->getFEEngine("SolidMechanicsFEEngine");
			Real average = 0;

			GhostType gt = _not_ghost;
			for (auto element_type :
			mesh.elementTypes(spatial_dimension, gt, _ek_not_defined)) {
			if (field_type == "stress") {
			for (UInt m = 0; m < this->materials.size(); ++m) {
			const auto & stress_vec = this->materials[m]->getStress(element_type);
			const auto & elem_filter =
			this->materials[m]->getElementFilter(element_type);
			Array<Real> int_stress_vec(elem_filter.size(),
			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.size(); ++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 auto & gradu_vec = this->materials[m]->getGradU(element_type);
			const auto & elem_filter =
			this->materials[m]->getElementFilter(element_type);
			Array<Real> int_gradu_vec(elem_filter.size(),
			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.size(); ++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 auto & eigen_gradu_vec =
			this->materials[m]->getInternal<Real>("eigen_grad_u")(element_type);
			const auto & elem_filter =
			this->materials[m]->getElementFilter(element_type);
			Array<Real> int_eigen_gradu_vec(elem_filter.size(),
			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.size(); ++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_ERROR("Averaging not implemented for this field!!!");
			}
			}

			return average / this->volume;

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::homogenizeStiffness(Matrix<Real> & C_macro) {
			AKANTU_DEBUG_IN();
			const Int 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)

			/// clear the eigenstrain
			Matrix<Real> zero_eigengradu(dim, dim, 0.);
			GhostType gt = _not_ghost;
			for (auto element_type : mesh.elementTypes(dim, gt, _ek_not_defined)) {
			auto & prestrain_vect =
			const_cast<Array<Real> &>(this->getMaterial("gel").getInternal<Real>(
			"eigen_grad_u")(element_type));
			auto prestrain_it =
			prestrain_vect.begin(spatial_dimension, spatial_dimension);
			auto prestrain_end =
			prestrain_vect.end(spatial_dimension, spatial_dimension);

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

			/// 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.);

			/// save the damage state before filling up cracks
			// ElementTypeMapReal saved_damage("saved_damage");
			// saved_damage.initialize(getFEEngine(), _nb_component = 1, _default_value =
			// 0);
			// this->fillCracks(saved_damage);

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

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

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

			/// drain cracks
			// this->drainCracks(saved_damage);
			/// compute effective stiffness
			Matrix<Real> eps_inverse(voigt_size, voigt_size);
			eps_inverse.inverse(strains);
			C_macro.mul<false, false>(stresses, eps_inverse);
			AKANTU_DEBUG_OUT();
			}

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

			auto & solver = this->getNonLinearSolver();
			solver.set("max_iterations", 2);
			solver.set("threshold", 1e-6);
			solver.set("convergence_type", SolveConvergenceCriteria::_solution);
			this->solveStep();

			/// 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");
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::homogenizeEigenGradU(
			Matrix<Real> & eigen_gradu_macro) {
			AKANTU_DEBUG_IN();
			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");
			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			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;

			auto tmp = std::make_shared<GelMaterialSelector>(*this, box_size, "gel",
			this->nb_gel_pockets, eps);
			tmp->setFallback(material_selector);
			material_selector = tmp;
			}

			this->assignMaterialToElements();
			// synchronize the element material arrays
			this->synchronize(SynchronizationTag::_constitutive_law_id);

			for (auto & material : materials) {
			/// init internals properties
			const auto type = material->getID();
			if (type.find("material_FE2") != std::string::npos)
			continue;
			material->initMaterial();
			}

			this->synchronize(SynchronizationTag::_clh_init_cl);

			if (this->non_local_manager) {
			this->non_local_manager->initialize();
			}
			// SolidMechanicsModel::initMaterials();

			AKANTU_DEBUG_OUT();
			}

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::fillCracks(ElementTypeMapReal & saved_damage) {
			const auto & mat_gel = this->getMaterial("gel");
			Real E_gel = mat_gel.get("E");
			Real E_homogenized = 0.;

			for (auto && mat : materials) {
			if (mat->getName() == "gel" \|\| mat->getName() == "FE2_mat")
			continue;

			Real E = mat->get("E");
			auto & damage = mat->getInternal<Real>("damage");

			for (auto && type : damage.elementTypes()) {
			const auto & elem_filter = mat->getElementFilter(type);
			auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type);

			auto sav_dam_it =
			make_view(saved_damage(type), nb_integration_point).begin();
			for (auto && data :
			zip(elem_filter, make_view(damage(type), nb_integration_point))) {
			auto el = std::get<0>(data);
			auto & dam = std::get<1>(data);
			Vector<Real> sav_dam = sav_dam_it[el];

			sav_dam = dam;

			for (auto q : arange(dam.size())) {
			E_homogenized = (E_gel - E) * dam(q) + E;
			dam(q) = 1. - (E_homogenized / E);
			}
			}
			}
			}
			}

			/* -------------------------------------------------------------------------- */
			void SolidMechanicsModelRVE::drainCracks(
			const ElementTypeMapReal & saved_damage) {
			for (auto && mat : materials) {
			if (mat->getName() == "gel" \|\| mat->getName() == "FE2_mat")
			continue;
			auto & damage = mat->getInternal<Real>("damage");

			for (auto && type : damage.elementTypes()) {
			const auto & elem_filter = mat->getElementFilter(type);
			auto nb_integration_point = getFEEngine().getNbIntegrationPoints(type);

			auto sav_dam_it =
			make_view(saved_damage(type), nb_integration_point).begin();
			for (auto && data :
			zip(elem_filter, make_view(damage(type), nb_integration_point))) {
			auto el = std::get<0>(data);
			auto & dam = std::get<1>(data);
			Vector<Real> sav_dam = sav_dam_it[el];

			dam = sav_dam;
			}
			}
			}
			}

			} // namespace akantu

rAKA/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc01d5a6530006bugfixes/65-Making-the…

solid_mechanics_model_RVE.ccView Options

rAKA/extra_packages/extra-materials/src/material_FE2/solid_mechanics_model_RVE.cc
01d5a6530006bugfixes/65-Making-the…

solid_mechanics_model_RVE.cc
View Options