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.org.
//
// 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

#include "../config/config.hpp"
#include "../linalg/linalg.hpp"
#include "fespace.hpp"
#include "gridfunc.hpp"
#include "linearform.hpp"
#include "bilininteg.hpp"
#include "staticcond.hpp"
#include "hybridization.hpp"

namespace mfem
{

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

   SparseMatrix *mat_e;

   FiniteElementSpace *fes;

   long sequence;

   int extern_bfs;

   Array<BilinearFormIntegrator*> dbfi;

   Array<BilinearFormIntegrator*> bbfi;

   Array<BilinearFormIntegrator*> fbfi;

   Array<BilinearFormIntegrator*> bfbfi;

   DenseMatrix elemmat;
   Array<int>  vdofs;

   DenseTensor *element_matrices;

   StaticCondensation *static_cond;
   Hybridization *hybridization;

   int precompute_sparsity;
   // Allocate appropriate SparseMatrix and assign it to mat
   void AllocMat();

   void ConformingAssemble();

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

public:
   BilinearForm(FiniteElementSpace *f);

   BilinearForm(FiniteElementSpace *f, BilinearForm *bf, int ps = 0);

   int Size() const { return height; }

   void EnableStaticCondensation();

   bool StaticCondensationIsEnabled() const { return static_cond; }

   FiniteElementSpace *SCFESpace() const
   { return static_cond ? static_cond->GetTraceFESpace() : NULL; }

   void EnableHybridization(FiniteElementSpace *constr_space,
                            BilinearFormIntegrator *constr_integ,
                            const Array<int> &ess_tdof_list);

   void UsePrecomputedSparsity(int ps = 1) { precompute_sparsity = ps; }

   void AllocateMatrix() { if (mat == NULL) { AllocMat(); } }

   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
   {
      MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
      return *mat;
   }
   SparseMatrix &SpMat()
   {
      MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
      return *mat;
   }
   SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }

   const SparseMatrix &SpMatElim() const
   {
      MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
      return *mat_e;
   }
   SparseMatrix &SpMatElim()
   {
      MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
      return *mat_e;
   }

   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 FormLinearSystem(Array<int> &ess_tdof_list, Vector &x, Vector &b,
                         SparseMatrix &A, Vector &X, Vector &B,
                         int copy_interior = 0);

   void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);

   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 AssembleBdrElementMatrix(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 EliminateEssentialBC(Array<int> &bdr_attr_is_ess, int d = 0);
   void EliminateEssentialBCDiag(Array<int> &bdr_attr_is_ess, double value);

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

   void EliminateVDofs(Array<int> &vdofs, 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 EliminateEssentialBCFromDofsDiag(Array<int> &ess_dofs, double value);

   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); }

   virtual 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;
   Array<BilinearFormIntegrator*> skt; // trace face integrators

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

   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; }
   SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }

   void AddDomainIntegrator (BilinearFormIntegrator * bfi);

   void AddBoundaryIntegrator (BilinearFormIntegrator * bfi);

   void AddTraceFaceIntegrator (BilinearFormIntegrator * bfi);

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

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

   Array<BilinearFormIntegrator*> *GetTFBFI() { return &skt; }

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

   void Assemble (int skip_zeros = 1);

   void ConformingAssemble();

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

   void EliminateEssentialBCFromTrialDofs(Array<int> &marked_vdofs,
                                          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); }

   void AddTraceFaceInterpolator(DiscreteInterpolator *di)
   { AddTraceFaceIntegrator(di); }

   Array<BilinearFormIntegrator*> *GetDI() { return &dom; }

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

}

#endif