bilinearform.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.googlecode.com.
//
// MFEM 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) version 2.1 dated February 1999.

#ifndef MFEM_BILINEARFORM
#define MFEM_BILINEARFORM


class BilinearForm : public Matrix
{
protected:
   SparseMatrix *mat;

   // Matrix used to eliminate b.c.
   SparseMatrix *mat_e;

   FiniteElementSpace *fes;

   int extern_bfs;

   Array<BilinearFormIntegrator*> dbfi;

   Array<BilinearFormIntegrator*> bbfi;

   Array<BilinearFormIntegrator*> fbfi;

   Array<BilinearFormIntegrator*> bfbfi;

   DenseMatrix elemmat;
   Array<int>  vdofs;

   DenseTensor *element_matrices;

   // may be used in the construction of derived classes
   BilinearForm() : Matrix (0)
   { fes = NULL; mat = mat_e = NULL; extern_bfs = 0; element_matrices = NULL; }

public:
   BilinearForm (FiniteElementSpace * f);

   BilinearForm (FiniteElementSpace * f, BilinearForm * bf);

   Array<BilinearFormIntegrator*> *GetDBFI() { return &dbfi; }

   Array<BilinearFormIntegrator*> *GetBBFI() { return &bbfi; }

   Array<BilinearFormIntegrator*> *GetFBFI() { return &fbfi; }

   Array<BilinearFormIntegrator*> *GetBFBFI() { return &bfbfi; }

   const double &operator()(int i, int j) { return (*mat)(i,j); }

   virtual double &Elem(int i, int j);

   virtual const double &Elem(int i, int j) const;

   virtual void Mult(const Vector &x, Vector &y) const;

   void FullMult(const Vector &x, Vector &y) const
   { mat->Mult(x, y); mat_e->AddMult(x, y); }

   virtual void AddMult(const Vector &x, Vector &y, const double a = 1.0) const
   { mat -> AddMult (x, y, a); }

   void FullAddMult(const Vector &x, Vector &y) const
   { mat->AddMult(x, y); mat_e->AddMult(x, y); }

   double InnerProduct(const Vector &x, const Vector &y) const
   { return mat->InnerProduct (x, y); }

   virtual MatrixInverse *Inverse() const;

   virtual void Finalize(int skip_zeros = 1);

   const SparseMatrix &SpMat() const { return *mat; }
   SparseMatrix &SpMat() { return *mat; }

   void AddDomainIntegrator(BilinearFormIntegrator *bfi);

   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi);

   void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi);

   void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi);

   void operator=(const double a)
   { if (mat != NULL) *mat = a; if (mat_e != NULL) *mat_e = a; }

   void Assemble(int skip_zeros = 1);

   void ComputeElementMatrices();

   void FreeElementMatrices()
   { delete element_matrices; element_matrices = NULL; }

   void ComputeElementMatrix(int i, DenseMatrix &elmat);
   void AssembleElementMatrix(int i, const DenseMatrix &elmat,
                              Array<int> &vdofs, int skip_zeros = 1);

   void EliminateEssentialBC(Array<int> &bdr_attr_is_ess,
                             Vector &sol, Vector &rhs, int d = 0);

   void EliminateVDofs(Array<int> &vdofs, Vector &sol, Vector &rhs, int d = 0);

   void EliminateVDofs(Array<int> &vdofs, int d = 0);

   void EliminateVDofsInRHS(Array<int> &vdofs, const Vector &x, Vector &b);

   double FullInnerProduct(const Vector &x, const Vector &y) const
   { return mat->InnerProduct(x, y) + mat_e->InnerProduct(x, y); }

   void EliminateEssentialBC(Array<int> &bdr_attr_is_ess, int d = 0);

   void EliminateEssentialBCFromDofs(Array<int> &ess_dofs, Vector &sol,
                                     Vector &rhs, int d = 0);

   void EliminateEssentialBCFromDofs(Array<int> &ess_dofs, int d = 0);

   void Update(FiniteElementSpace *nfes = NULL);

   FiniteElementSpace *GetFES() { return fes; }

   virtual ~BilinearForm();
};

class MixedBilinearForm : public Matrix
{
protected:
   SparseMatrix *mat;

   FiniteElementSpace *trial_fes, *test_fes;

   Array<BilinearFormIntegrator*> dom;
   Array<BilinearFormIntegrator*> bdr;

   int width;

public:
   MixedBilinearForm (FiniteElementSpace *tr_fes,
                      FiniteElementSpace *te_fes);

   int Width() const { return width; };

   virtual double& Elem (int i, int j);

   virtual const double& Elem (int i, int j) const;

   virtual void Mult (const Vector & x, Vector & y) const;

   virtual void AddMult (const Vector & x, Vector & y,
                         const double a = 1.0) const;

   virtual void AddMultTranspose (const Vector & x, Vector & y,
                                  const double a = 1.0) const;

   virtual void MultTranspose (const Vector & x, Vector & y) const
   { y = 0.0; AddMultTranspose (x, y); };

   virtual MatrixInverse * Inverse() const;

   virtual void Finalize (int skip_zeros = 1);

   void GetBlocks(Array2D<SparseMatrix *> &blocks) const;

   const SparseMatrix &SpMat() const { return *mat; };
   SparseMatrix &SpMat() { return *mat; };

   void AddDomainIntegrator (BilinearFormIntegrator * bfi);

   void AddBoundaryIntegrator (BilinearFormIntegrator * bfi);

   void operator= (const double a) { *mat = a; }

   virtual void Assemble (int skip_zeros = 1);

   void EliminateTrialDofs (Array<int> &bdr_attr_is_ess,
                            Vector &sol, Vector &rhs);

   virtual void EliminateTestDofs (Array<int> &bdr_attr_is_ess);

   void Update();

   virtual ~MixedBilinearForm();
};


class DiscreteLinearOperator : public MixedBilinearForm
{
public:
   DiscreteLinearOperator(FiniteElementSpace *domain_fes,
                          FiniteElementSpace *range_fes)
      : MixedBilinearForm(domain_fes, range_fes) { }

   void AddDomainInterpolator(DiscreteInterpolator *di)
   { AddDomainIntegrator(di); }

   virtual void Assemble(int skip_zeros = 1);
};

#endif