Added Newton nonlinear solver. Still does not work.

AMDiS/src/nonlin/Newton.h

+// ============================================================================
+// ==                                                                        ==
+// == AMDiS - Adaptive multidimensional simulations                          ==
+// ==                                                                        ==
+// ==  http://www.amdis-fem.org                                              ==
+// ==                                                                        ==
+// ============================================================================
+//
+// Software License for AMDiS
+//
+// Copyright (c) 2010 Dresden University of Technology 
+// All rights reserved.
+// Authors: Simon Vey, Thomas Witkowski et al.
+//
+// This file is part of AMDiS
+//
+// See also license.opensource.txt in the distribution.
+
+
+
+/** \file Newton.h */
+
+#ifndef AMDIS_NEWTON_H
+#define AMDIS_NEWTON_H
+
+#include "CreatorInterface.h"
+#include "NonLinSolver.h"
+#include "OEMSolver.h"
+
+namespace AMDiS {
+
+  /**
+   * \ingroup Solver
+   * 
+   * \Brief
+   * Implements the newton method for solving a non linear system. Sub class of
+   * NonLinSolver.
+   */
+  class Newton : public NonLinSolver
+  {
+  public:
+    /// Creator class used in the NonLinSolverMap.
+    class Creator : public NonLinSolverCreator
+    {
+    public:
+      virtual ~Creator() {}
+
+      /** \brief
+       * Returns a new Newton object.
+       */
+      NonLinSolver* create() 
+      { 
+	return new Newton(this->name, this->linearSolver, this->nonLinUpdater); 
+      }
+    };
+
+    /** \brief
+     * Calls constructor of base class NonLinSolver
+     */
+    Newton(const std::string& name_, 
+	   OEMSolver *linSolver_,
+	   NonLinUpdater *updater)
+      : NonLinSolver(name_, linSolver_, updater),
+	b(NULL)
+    {}
+
+  private:
+    /** \brief
+     * realisation of NonLinSolver::init
+     */
+    void init()
+    {}
+
+    /** \brief
+     * realisation of NonLinSolver::nlsolve
+     */
+    int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& mat,
+		SystemVector& x, SystemVector& rhs)
+    {
+      //DOFVector<T> *d = problem->getRHS();
+      //DOFVector<T> *x = problem->getSolution();;
+
+      b = new SystemVector(x);
+      *b = rhs;
+
+      //   // copy operators from fh to b
+      //   std::vector<Operator*>::iterator op;
+      //   std::vector<double*>::iterator fac;
+      //   for(op = d->getOperatorsBegin(), 
+      // 	fac = d->getOperatorFactorBegin();
+      //       op != d->getOperatorsEnd();
+      //       ++op, ++fac) 
+      //   {
+      //     b->addOperator(*op, *fac);
+      //   }
+
+      double err = 0.0, err_old = -1.0;
+      int    iter, n;
+
+      INFO(this->info,2)("iter. |     this->residual |     red. |    n |\n");
+
+      for (iter = 1; iter <= this->maxIter; iter++) {
+	/*--- Assemble DF(x) and F(x) ----------------------------------------------*/
+	this->nonLinUpdater->update(/*x,*/true, b);
+	/*--- Initial guess is zero ------------------------------------------------*/
+	rhs.set(0.0);
+	/*--- solve linear system --------------------------------------------------*/
+	n = solveLinearSystem(mat, *b, rhs);
+	/*--- x = x - d ------------------------------------------------------------*/
+	x -= rhs;
+
+	if (this->usedNorm == NO_NORM || this->usedNorm == L2_NORM)
+	  err = L2Norm(&rhs);  // sollte hier nicht b genommen werden (s. NewtonS) ?
+	else
+	  err = H1Norm(&rhs);  // sollte hier nicht b genommen werden (s. NewtonS) ?
+    
+
+	if (iter == 1)  this->initial_residual = err;
+
+	if (err_old <= 0) {
+	  INFO(this->info,2)("%5d | %12.5e | -------- | %4d |\n", iter, err, n);
+	} else {
+	  INFO(this->info,2)("%5d | %12.5e | %8.2e | %4d |\n", 
+			     iter, err, err/err_old, n);
+	}
+
+	if ((this->residual = err) < this->tolerance) {
+	  INFO(this->info,4)("finished successfully\n");
+	  return iter;
+	}
+	err_old = err;
+      }
+
+      if (this->info < 2)
+	INFO(this->info,1)("iter. %d, residual: %12.5e\n", iter, err);
+      INFO(this->info,1)("tolerance %e not reached\n", this->tolerance);
+
+      this->residual = err;
+
+      return iter;
+    }
+
+    /** \brief
+     * realisation of NonLinSolver::exit
+     */
+    void exit()
+    {
+      if (b != NULL)
+	delete b;
+    }
+
+  private:
+    /** \brief
+     * internal used data
+     */
+    SystemVector *b;
+  };
+
+}
+
+#endif // AMDIS_NEWTON_H

AMDiS/src/nonlin/NonLinSolver.h

@@ -181,5 +181,7 @@ namespace AMDiS {
 
 }
 
+#include "Newton.h"
+
 #endif // AMDIS_NONLINSOLVER_H