Homec4science
Diffusion akantu (features/fe2-orthodamage) Browse

rAKA/src/model/multiscale_tools/solid_mechanics_model_RVE.cc
2d7ec870588efeatures/fe2-orthodamage

Last Change
Raw File

solid_mechanics_model_RVE.cc
View Options

/**
* @file solid_mechanics_model_RVE.cc
* @author Emil Gallyamov <emil.gallyamov@epfl.ch>
* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Wed Jan 13 15:32:35 2016
* @update Tue Feb 8 2022
*
* @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 "aka_random_generator.hh"
#include "element_group.hh"
#include "expanding_material_selector.hh"
#include "material_damage_iterative.hh"
// #include "model_options.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, UInt nb_expanding_elements, bool use_RVE_mat_selector,
UInt dim, const ID & id, std::shared_ptr<DOFManager> dof_manager)
: SolidMechanicsModel(mesh, dim, id, dof_manager,
ModelType::_solid_mechanics_model_rve),
RVETools(dynamic_cast<SolidMechanicsModel &>(*this)),
use_RVE_mat_selector(use_RVE_mat_selector),
nb_expanding_elements(nb_expanding_elements), stiffness_changed(true) {
AKANTU_DEBUG_IN();
// seed is fixed to the same value for all the processors
// It is done to generate same random internal fields
RandomGenerator<UInt>::seed(2);
// 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());
mesh.makePeriodic(_x, "left", "right");
mesh.makePeriodic(_y, "bottom", "top");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModelRVE::~SolidMechanicsModelRVE() = default;
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::initFullImpl(const ModelOptions & options) {
AKANTU_DEBUG_IN();
SolidMechanicsModel::initFullImpl(options);
auto & solver = this->getNonLinearSolver();
solver.set("max_iterations", 50);
solver.set("threshold", 1e-4);
solver.set("convergence_type", SolveConvergenceCriteria::_solution);
// compute the volume of the RVE
this->computeModelVolume();
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->addDumpField("extra_volume");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelRVE::assembleInternalForces() {
// displacement correction
auto disp_begin = make_view(*this->displacement, spatial_dimension).begin();
auto u_top_corr = Vector<Real>(disp_begin[this->corner_nodes(2)]) -
Vector<Real>(disp_begin[this->corner_nodes(1)]);
auto u_right_corr = Vector<Real>(disp_begin[this->corner_nodes(2)]) -
Vector<Real>(disp_begin[this->corner_nodes(3)]);
auto correct_disp = [&](auto && node, auto && correction) {
const auto & pair = mesh.getPeriodicMasterSlaves().equal_range(node);
for (auto && data : range(pair.first, pair.second)) {
auto slave = data.second;
auto slave_disp = Vector<Real>(disp_begin[slave]);
slave_disp = Vector<Real>(disp_begin[node]);
slave_disp += correction;
}
};
for (auto node : bottom_nodes) {
correct_disp(node, u_top_corr);
}
for (auto node : left_nodes) {
correct_disp(node, u_right_corr);
}
SolidMechanicsModel::assembleInternalForces();
}
/* --------------------------------------------------------------------------
*/
void SolidMechanicsModelRVE::advanceExpansion(const Matrix<Real> & prestrain,
const ID & material_name) {
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(material_name)
.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
MaterialDamageIterativeInterface & mat_paste =
dynamic_cast<MaterialDamageIterativeInterface &>(
this->getMaterial("paste"));
MaterialDamageIterativeInterface & mat_aggregate =
dynamic_cast<MaterialDamageIterativeInterface &>(
this->getMaterial("aggregate"));
UInt nb_damaged_elements = 0;
Real max_eq_stress_aggregate = 0;
Real max_eq_stress_paste = 0;
this->stiffness_changed = false;
// variables for parallel execution
auto && comm = akantu::Communicator::getWorldCommunicator();
auto prank = comm.whoAmI();
// save nodal fields (disp, boun & ext_force)
this->storeNodalFields();
do {
// restore nodals and update grad_u accordingly
this->restoreNodalFields();
// restore historical internal fields (sigma_v for visc)
this->restoreInternalFields();
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());
}
// mark the flag to update stiffness if elements were damaged
if (nb_damaged_elements != 0) {
this->stiffness_changed = true;
}
AKANTU_DEBUG_INFO("Proc " << prank << " the number of damaged elements is "
<< nb_damaged_elements);
} while (nb_damaged_elements);
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) {
auto tmp = std::make_shared<ExpandingMaterialSelector>(
*this, "gel", this->nb_expanding_elements, "aggregate");
tmp->setFallback(material_selector);
material_selector = tmp;
}
this->assignMaterialToElements();
// synchronize the element material arrays
this->synchronize(SynchronizationTag::_material_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::_smm_init_mat);
if (this->non_local_manager) {
this->non_local_manager->initialize();
}
AKANTU_DEBUG_OUT();
}
} // namespace akantu
c4science · Help