Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/continuum/akantu/domain_akantu_static.cc
b05d5f2aba78features/reshape_comput…

Last Change
Raw File

domain_akantu_static.cc
View Options

/**
* @file domain_akantu_static.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Jul 23 19:06:50 2014
*
* @brief This is the model wrapping Akantu in its static solve version
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#define TIMER
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#include "domain_akantu_static.hh"
#include "lib_geometry.hh"
#include "math_tools.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <aka_types.hh>
#include <solid_mechanics_model.hh>
#include <material_elastic.hh>
/* -------------------------------------------------------------------------- */
#include <iomanip>
#include <sstream>
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <UInt Dim>
DomainAkantuStatic<Dim>::~DomainAkantuStatic(){
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
DomainAkantuStatic<Dim>::DomainAkantuStatic(DomainID ID,CommGroup GID):
DomainAkantu<Dim>(ID,GID), first_solve(true)
{}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuStatic<Dim>::init(){
if (this->material_filename != "")
this->akantu_parser.parse(this->material_filename);
this->model->setParser(this->akantu_parser);
this->model->initFull(akantu::SolidMechanicsModelOptions(akantu::_static,true));
this->initMaterial();
this->elems.setSolidMechanicsModel(*this->model);
this->nodes.setSolidMechanicsModel(*this->model);
this->model->setTimeStep(-1);
if (!this->nonlinear) {
this->model->assembleStiffnessMatrix();
}
this->model->assembleMassLumped();
akantu::Array<Real>::const_vector_iterator mass_it = this->model->getMass().begin(Dim);
akantu::Array<Real>::const_vector_iterator mass_end = this->model->getMass().end(Dim);
akantu::Vector<Real> sum(Dim, 0.);
UInt nb_nodes = this->model->getMesh().getNbNodes();
for (UInt n = 0 ; n < nb_nodes ; ++n, ++mass_it) {
if(!this->model->isPBCSlaveNode(n) && this->model->getMesh().isLocalOrMasterNode(n)) {
sum += *mass_it;
}
}
this->mass_sum = 0;
for (UInt i = 0 ; i < Dim ; ++i) {
this->mass_sum += sum(i)/Dim;
}
this->mass_sum = sum.norm<akantu::L_1>() / Dim;
DUMP("Total mass of mesh: " << this->mass_sum, DBG_MESSAGE);
//this->model->addDumpFieldVector("displacement");
//this->model->addDumpFieldVector("residual");
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuStatic<Dim>::performStep1(){
STARTTIMER("AkantuStepPosition");
DomainAkantu<Dim>::performStep1();
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuStatic<Dim>::performStep2(){
Real error = 0, initial_error = 0;
UInt niter = 0;
if (!LMsolveStep(this->tolerance, this->max_iter, error, initial_error, niter)) {
LM_FATAL(this->getID() << ": Static Newton-Raphson of akantu did not converge with a tolerance of " << this->tolerance << " per unit mass. It failed with an error norm of " << error << " per unit mass (initially " << initial_error << ") after " << niter << " iterations.");
} else if (niter > 0 ) {
DUMP(this->getID() << ": Static Newton-Raphson converged after " << niter << " iterations with an error norm of " << error << " per unit mass (initially " << initial_error << "). (Tolerance = " << this->tolerance << " per unit mass).", DBG_MESSAGE);//INFO_STARTUP);
} else {
DUMP(this->getID() << ": Static Newton-Raphson converged after " << niter << " iterations with an error norm of " << error << " per unit mass (initially " << initial_error << "). (Tolerance = " << this->tolerance << " per unit mass).", DBG_DETAIL);
}
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuStatic<Dim>::LMupdateResidual() {
this->model->initializeUpdateResidualData();
if (this->prestressed_periodicity_flag) {
typedef std::vector<std::pair<UInt, UInt> >::iterator pair_it;
pair_it pbc_pair = this->pbc_pairs.begin();
pair_it end = this->pbc_pairs.end();
UInt node_counter = 0;
for (; pbc_pair != end; ++pbc_pair, ++node_counter) {
for (UInt direction = 0 ; direction < Dim ; ++direction) {
(*this->displacement)[Dim*pbc_pair->first + direction] +=
this->slave_displacements[Dim * node_counter+direction];
}
}
}
this->model->updateResidual(false);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
bool DomainAkantuStatic<Dim>::LMsolveStep(Real tol, UInt max_iteration, Real & error,
Real & initial_error, UInt & niter) {
if (this->first_solve) {
akantu::Array<bool> & bound = this->model->getBlockedDOFs();
UInt full_counter = 0;
for (UInt i = 0 ; i < bound.getSize() ; ++i) {
for (UInt j = 0 ; j < Dim ; ++j) {
if (bound(i, j)) { // yeah, yeah, I know...
++full_counter;
}
}
}
UInt min_nb_bound = 1;
if (Dim == 2) min_nb_bound = 3;
else if (Dim == 3) min_nb_bound = 6;
Communicator::getCommunicator().allReduce(&full_counter,1,this->getGroupID(),"full_counter",OP_SUM);
if (full_counter < min_nb_bound) {
LM_FATAL("The number of blocked dofs is insufficient. You specified "
<< full_counter << " constraints, and a " << Dim << "D problem needs "
<< "at least " << min_nb_bound << " blocked dofs!");
}
}
this->model->implicitPred();
this->LMupdateResidual();
bool converged = false;
if (!this->first_solve) {
converged = this->model->template testConvergence<akantu::_scc_residual_mass_wgh> (tol, initial_error);
error = initial_error;
if(converged) {
return converged;
}
}
do {
if (this->nonlinear) {
this->model->assembleStiffnessMatrix();
}
this->model->template solve<akantu::NewmarkBeta::_displacement_corrector> (this->model->getIncrement(), 1.,
this->first_solve||this->nonlinear);
this->first_solve = false;
this->model->implicitCorr();
//this->model->dump();
this->LMupdateResidual();
//this->model->dump();
converged = this->model-> template testConvergence<akantu::_scc_residual_mass_wgh> (tol, error);
niter++;
DUMP("[" << akantu::_scc_residual << "] Convergence iteration "
<< std::setw(std::log10(max_iteration) + 1) << niter
<< ": error " << error << (converged ? " < " : " > ") << tol,
DBG_DETAIL);
} while (!converged && niter < max_iteration);
if (!converged || current_step ==1) {
this->model->getGlobalJacobianMatrix().saveMatrix("matrix_at_failure.mtx");
}
return converged;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void DomainAkantuStatic<Dim>::performStep3(){
this->mesh_container.incRelease();
}
/* -------------------------------------------------------------------------- */
/* LMDESC AKANTU_STATIC
This domain implements the plugin with Akantu Finite Element library.
*/
/* LMHERITANCE domain_akantu */
template <UInt Dim>
void DomainAkantuStatic<Dim>::declareParams(){
//DomainAkantu<Dim>::declareParams();
/* LMKEYWORD NONLINEAR
specifies whether the continuum material is nonlinear, in which case
the stiffness matrix is reassembled at every Newton-Raphson iteration
*/
this->parseKeyword("NONLINEAR", this->nonlinear, false);
/* LMKEYWORD TOLERANCE
residual tolerance for Newton-Raphson loop. This tolerance is multiplied by
total mass of the mesh to make it scale. So think in force-tolerance per
mass when choosing this value
*/
this->parseKeyword("TOLERANCE", this->tolerance, 1e-3);
/* LMKEYWORD MAX_ITER
max iterations for Newton-Raphson loop. for linear materials 1 should
be enough
*/
this->parseKeyword("MAX_ITER", this->max_iter, 1u);
}
/* -------------------------------------------------------------------------- */
template class DomainAkantuStatic<1>;
template class DomainAkantuStatic<2>;
template class DomainAkantuStatic<3>;
__END_LIBMULTISCALE__
c4science · Help