bilininteg.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_BILININTEG
#define MFEM_BILININTEG

class BilinearFormIntegrator
{
public:
   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
   virtual void AssembleFaceMatrix(const FiniteElement &el1,
                                   const FiniteElement &el2,
                                   FaceElementTransformations &Trans,
                                   DenseMatrix &elmat);
   virtual void ComputeElementFlux(const FiniteElement &el,
                                   ElementTransformation &Trans,
                                   Vector &u,
                                   const FiniteElement &fluxelem,
                                   Vector &flux, int wcoef = 1) { }
   virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                    ElementTransformation &Trans,
                                    Vector &flux) { return 0.0; }
   virtual ~BilinearFormIntegrator() { }
};

class TransposeIntegrator : public BilinearFormIntegrator
{
private:
   int own_bfi;
   BilinearFormIntegrator *bfi;

   DenseMatrix bfi_elmat;

public:
   TransposeIntegrator (BilinearFormIntegrator *_bfi, int _own_bfi = 1)
   { bfi = _bfi; own_bfi = _own_bfi; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
   virtual ~TransposeIntegrator() { if (own_bfi) delete bfi; }
};

class LumpedIntegrator : public BilinearFormIntegrator
{
private:
   int own_bfi;
   BilinearFormIntegrator *bfi;

public:
   LumpedIntegrator (BilinearFormIntegrator *_bfi, int _own_bfi = 1)
   { bfi = _bfi; own_bfi = _own_bfi; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);

   virtual ~LumpedIntegrator() { if (own_bfi) delete bfi; }
};

class DiffusionIntegrator: public BilinearFormIntegrator
{
private:
   Vector vec, pointflux, shape;
#ifndef MFEM_USE_OPENMP
   DenseMatrix dshape, dshapedxt, invdfdx;
#endif
   Coefficient *Q;
   MatrixCoefficient *MQ;

public:
   DiffusionIntegrator (Coefficient &q) : Q(&q) { MQ = NULL; }

   DiffusionIntegrator (MatrixCoefficient &q) : MQ(&q) { Q = NULL; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
   virtual void ComputeElementFlux(const FiniteElement &el,
                                   ElementTransformation &Trans,
                                   Vector &u,
                                   const FiniteElement &fluxelem,
                                   Vector &flux, int wcoef);
   virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                    ElementTransformation &Trans,
                                    Vector &flux);
};

class MassIntegrator: public BilinearFormIntegrator
{
private:
   Vector shape, te_shape;
   Coefficient *Q;
   const IntegrationRule *IntRule;

public:
   MassIntegrator(const IntegrationRule *ir = NULL)
   { Q = NULL; IntRule = ir; }
   MassIntegrator(Coefficient &q, const IntegrationRule *ir = NULL)
      : Q(&q) { IntRule = ir; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
};

class BoundaryMassIntegrator : public MassIntegrator
{
public:
   BoundaryMassIntegrator(Coefficient &q) : MassIntegrator(q) { }
};

class ConvectionIntegrator : public BilinearFormIntegrator
{
private:
   DenseMatrix dshape;
   DenseMatrix invdfdx;
   Vector shape, vec1, vec2;
   Vector BdFidxT;
   VectorCoefficient &Q;
public:
   ConvectionIntegrator(VectorCoefficient &q) : Q(q) { }
   virtual void AssembleElementMatrix(const FiniteElement &,
                                      ElementTransformation &,
                                      DenseMatrix &);
};

class VectorMassIntegrator: public BilinearFormIntegrator
{
private:
   Vector shape, vec;
   DenseMatrix partelmat;
   DenseMatrix mcoeff;
   Coefficient *Q;
   VectorCoefficient *VQ;
   MatrixCoefficient *MQ;

   const IntegrationRule *IntRule;
   int Q_order;

public:
   VectorMassIntegrator()
   { Q = NULL; VQ = NULL; MQ = NULL; IntRule = NULL; Q_order = 0; }
   VectorMassIntegrator(Coefficient &q, int qo = 0)
      : Q(&q) { VQ = NULL; MQ = NULL; IntRule = NULL; Q_order = qo; }
   VectorMassIntegrator(Coefficient &q, const IntegrationRule *ir)
      : Q(&q) { VQ = NULL; MQ = NULL; IntRule = ir; Q_order = 0; }
   VectorMassIntegrator(VectorCoefficient &q, int qo = 0)
      : VQ(&q) { Q = NULL; MQ = NULL; IntRule = NULL; Q_order = qo; }
   VectorMassIntegrator(MatrixCoefficient &q, int qo = 0)
      : MQ(&q) { Q = NULL; VQ = NULL; IntRule = NULL; Q_order = qo; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
};


class VectorFEDivergenceIntegrator : public BilinearFormIntegrator
{
private:
   Coefficient *Q;
#ifndef MFEM_USE_OPENMP
   Vector divshape, shape;
#endif
public:
   VectorFEDivergenceIntegrator() { Q = NULL; }
   VectorFEDivergenceIntegrator(Coefficient &q) { Q = &q; }
   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat) { }
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
};


class VectorFECurlIntegrator: public BilinearFormIntegrator
{
private:
   Coefficient *Q;
#ifndef MFEM_USE_OPENMP
   DenseMatrix curlshapeTrial;
   DenseMatrix vshapeTest;
   DenseMatrix curlshapeTrial_dFT;
#endif
public:
   VectorFECurlIntegrator() { Q = NULL; }
   VectorFECurlIntegrator(Coefficient &q) { Q = &q; }
   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat) { }
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
};


class DerivativeIntegrator : public BilinearFormIntegrator
{
private:
   Coefficient & Q;
   int xi;
   DenseMatrix dshape, dshapedxt, invdfdx;
   Vector shape, dshapedxi;
public:
   DerivativeIntegrator(Coefficient &q, int i) : Q(q), xi(i) { }
   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat)
   { AssembleElementMatrix2(el,el,Trans,elmat); }
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
};

class CurlCurlIntegrator: public BilinearFormIntegrator
{
private:
#ifndef MFEM_USE_OPENMP
   DenseMatrix Curlshape, Curlshape_dFt;
#endif
   Coefficient *Q;

public:
   CurlCurlIntegrator() { Q = NULL; }
   CurlCurlIntegrator(Coefficient &q) : Q(&q) { }

   /* Given a particular Finite Element, compute the
      element curl-curl matrix elmat */
   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
};

class VectorFEMassIntegrator: public BilinearFormIntegrator
{
private:
   Coefficient *Q;
   VectorCoefficient *VQ;

#ifndef MFEM_USE_OPENMP
   Vector shape;
   Vector D;
   DenseMatrix vshape;
#endif

public:
   VectorFEMassIntegrator() { Q = NULL; VQ= NULL; }
   VectorFEMassIntegrator(Coefficient *_q) { Q = _q; VQ= NULL; }
   VectorFEMassIntegrator(Coefficient &q) { Q = &q; VQ= NULL; }
   VectorFEMassIntegrator(VectorCoefficient *_vq) { VQ = _vq; Q = NULL; }
   VectorFEMassIntegrator(VectorCoefficient &vq) { VQ = &vq; Q = NULL; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
};

class VectorDivergenceIntegrator : public BilinearFormIntegrator
{
private:
   Coefficient *Q;

   Vector shape;
   Vector divshape;
   DenseMatrix dshape;
   DenseMatrix gshape;
   DenseMatrix Jadj;

public:
   VectorDivergenceIntegrator() { Q = NULL; }
   VectorDivergenceIntegrator(Coefficient *_q) { Q = _q; }
   VectorDivergenceIntegrator(Coefficient &q) { Q = &q; }

   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
                                       const FiniteElement &test_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
};

class DivDivIntegrator: public BilinearFormIntegrator
{
private:
   Coefficient *Q;

#ifndef MFEM_USE_OPENMP
   Vector divshape;
#endif

public:
   DivDivIntegrator() { Q = NULL; }
   DivDivIntegrator(Coefficient &q) : Q(&q) { }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
};

class VectorDiffusionIntegrator : public BilinearFormIntegrator
{
private:
   Coefficient *Q;

   DenseMatrix Jinv;
   DenseMatrix dshape;
   DenseMatrix gshape;
   DenseMatrix pelmat;

public:
   VectorDiffusionIntegrator() { Q = NULL; }
   VectorDiffusionIntegrator(Coefficient &q) { Q = &q; }

   virtual void AssembleElementMatrix(const FiniteElement &el,
                                      ElementTransformation &Trans,
                                      DenseMatrix &elmat);
};

class ElasticityIntegrator : public BilinearFormIntegrator
{
private:
   double q_lambda, q_mu;
   Coefficient *lambda, *mu;

#ifndef MFEM_USE_OPENMP
   DenseMatrix dshape, Jinv, gshape, pelmat;
   Vector divshape;
#endif

public:
   ElasticityIntegrator(Coefficient &l, Coefficient &m)
   { lambda = &l; mu = &m; }
   ElasticityIntegrator(Coefficient &m, double q_l, double q_m)
   { lambda = NULL; mu = &m; q_lambda = q_l; q_mu = q_m; }

   virtual void AssembleElementMatrix(const FiniteElement &,
                                      ElementTransformation &,
                                      DenseMatrix &);
};


class DiscreteInterpolator : public BilinearFormIntegrator { };


class GradientInterpolator : public DiscreteInterpolator
{
public:
   virtual void AssembleElementMatrix2(const FiniteElement &h1_fe,
                                       const FiniteElement &nd_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat)
   { nd_fe.ProjectGrad(h1_fe, Trans, elmat); }
};


class IdentityInterpolator : public DiscreteInterpolator
{
public:
   virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
                                       const FiniteElement &ran_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat)
   { ran_fe.Project(dom_fe, Trans, elmat); }
};


class CurlInterpolator : public DiscreteInterpolator
{
public:
   virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
                                       const FiniteElement &ran_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat)
   { ran_fe.ProjectCurl(dom_fe, Trans, elmat); }
};


class DivergenceInterpolator : public DiscreteInterpolator
{
public:
   virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
                                       const FiniteElement &ran_fe,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat)
   { ran_fe.ProjectDiv(dom_fe, Trans, elmat); }
};

#endif