Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/continuum/akantu/domain_akantu_static.cc
1e97b3b33259master

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 "domain_akantu_static.hh"
#include "lib_geometry.hh"
#include "lm_common.hh"
#include "math_tools.hh"
#include "units.hh"
/* -------------------------------------------------------------------------- */
#include <aka_types.hh>
#include <material_elastic.hh>
#include <solid_mechanics_model.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