Homec4science
Diffusion akantu (bugfix/69-segmentation-fault-caused-by-model-applybc-aka-fromtraction-trac-traction) Browse

rAKA/src/material_FE2/solid_mechanics_model_RVE.cc
0cfe38a37116bugfix/69-segmentation…

Last Change
Raw File

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__
c4science · Help