NewtonIteration.inc.hpp

#pragma once
 
#include <amdis/Integrate.hpp>
#include <amdis/operations/Arithmetic.hpp>
 
namespace AMDiS {
 
template <class Problem>
Flag NewtonIteration<Problem>::oneIteration(AdaptInfo& adaptInfo, Flag toDo)
{
  if (!toDo.isSet(SOLVE))
    return Flag(1);
 
  int iter = adaptInfo.spaceIteration()+1;
  if (iter <= 0 || (buildCycle_ > 0 && (iter-1) % buildCycle_ == 0)) {
    prob_->buildAfterAdapt(adaptInfo, 0, true, true);  // assemble DF(u0) and F(u)
    prob_->solver()->init(prob_->systemMatrix()->impl());
  }
  else {
    prob_->buildAfterAdapt(adaptInfo, 0, false, true); // assemble F(u)
  }
 
  // Initial guess is zero
  stepSolution_->resizeZero();
 
  // solve the linear system
  Dune::InverseOperatorResult stat;
  prob_->solver()->apply(stepSolution_->impl(), prob_->rhsVector()->impl(), stat);
 
  // u = u + d
  Recursive::transform(prob_->solutionVector()->impl(),
    Operation::Plus{}, prob_->solutionVector()->impl(), stepSolution_->impl());
 
  // L2-norm
  err_ = integrate(sqr(this->stepSolution()), prob_->gridView());
  if (norm_ > 1) // H1-norm
    err_ += integrate(unary_dot(gradientOf(this->stepSolution())), prob_->gridView());
  err_ = std::sqrt(err_);
 
  adaptInfo.setEstSum(err_, "");
  return Flag(1);
}
 
} // end namespace AMDiS