ex0_nonlinear.cpp

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

static const double alpha = -0.5;

double f(double u) { return exp(alpha*u); }
double df(double u) { return alpha*exp(alpha*u); }

// Define a coefficient that, given a grid function u, returns f(u)
class NonlinearCoefficient : public Coefficient
{
   GridFunction &gf;
public:
   NonlinearCoefficient(GridFunction &gf_) : gf(gf_) { }
   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
   {
      return f(gf.GetValue(T, ip));
   }
};

// Define a coefficient that, given a grid function u, returns df(u)
class NonlinearDerivativeCoefficient : public Coefficient
{
   GridFunction &gf;
public:
   NonlinearDerivativeCoefficient(GridFunction &gf_) : gf(gf_) { }
   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
   {
      return df(gf.GetValue(T, ip));
   }
};

// Define a nonlinear integrator that computes (f(u), v) and its linearized
// operator, (u df(u), v).
//
// Note that the action (f(u), v) can be computed using DomainLFIntegrator
// and the Jacobian matrix linearized operator can be computed using
// MassIntegrator with the appropriate coefficients.
class NonlinearMassIntegrator : public NonlinearFormIntegrator
{
   FiniteElementSpace &fes;
   GridFunction gf;
   Array<int> dofs;

public:
   NonlinearMassIntegrator(FiniteElementSpace &fes_) : fes(fes_), gf(&fes) { }

   virtual void AssembleElementVector(const FiniteElement &el,
                                      ElementTransformation &Tr,
                                      const Vector &elfun, Vector &elvect)
   {
      fes.GetElementDofs(Tr.ElementNo, dofs);
      gf.SetSubVector(dofs, elfun);
      NonlinearCoefficient coeff(gf);
      DomainLFIntegrator integ(coeff);
      Vector elvect_contrib;
      integ.AssembleRHSElementVect(el, Tr, elvect);
   }

   virtual void AssembleElementGrad(const FiniteElement &el,
                                    ElementTransformation &Tr,
                                    const Vector &elfun, DenseMatrix &elmat)
   {
      fes.GetElementDofs(Tr.ElementNo, dofs);
      gf.SetSubVector(dofs, elfun);
      NonlinearDerivativeCoefficient coeff(gf);
      MassIntegrator integ(coeff);
      integ.AssembleElementMatrix(el, Tr, elmat);
   }
};

int main(int argc, char *argv[])
{
   // 1. Parse command line options
   const char *mesh_file = "../data/star.mesh";
   int order = 1;

   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
   args.ParseCheck();

   // 2. Read the mesh from the given mesh file, and refine once uniformly.
   Mesh mesh(mesh_file);
   mesh.UniformRefinement();

   // 3. Define a finite element space on the mesh. Here we use H1 continuous
   //    high-order Lagrange finite elements of the given order.
   H1_FECollection fec(order, mesh.Dimension());
   FiniteElementSpace fespace(&mesh, &fec);
   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;

   // 4. Extract the list of all the boundary DOFs. These will be marked as
   //    Dirichlet in order to enforce zero boundary conditions.
   Array<int> boundary_dofs;
   fespace.GetBoundaryTrueDofs(boundary_dofs);

   // 5. Define the solution x as a finite element grid function in fespace. Set
   //    the initial guess to zero, which also sets the boundary conditions.
   GridFunction x(&fespace);
   x = 0.0;

   // 6. Set up the nonlinear form n(u,v) = (grad u, grad v) + (f(u), v)
   NonlinearForm n(&fespace);
   n.AddDomainIntegrator(new NonlinearMassIntegrator(fespace));
   n.AddDomainIntegrator(new DiffusionIntegrator);
   n.SetEssentialBC(boundary_dofs);

   // 7. Set up the the right-hand side. For simplicitly, we just use a zero
   //    vector. Because of the form of the nonlinear function f, it is still
   //    nontrivial to solve n(u,v) = 0.
   LinearForm b(&fespace);
   b = 0.0;

   // 8. Set up the Newton solver. Each Newton iteration requires a linear
   //    solve. Here we use UMFPack as a direct solver for these systems.
   UMFPackSolver direct_solver;
   NewtonSolver newton;
   newton.SetOperator(n);
   newton.SetSolver(direct_solver);
   newton.SetPrintLevel(1);
   newton.SetRelTol(1e-10);
   newton.SetMaxIter(20);

   // 9. Solve the nonlinear system.
   newton.Mult(b, x);

   return 0;
}