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