app-bholger/gfem/gfemLaplace.h

Go to the documentation of this file.

#pragma once

#include "gfemQuad.h"
#include "gfemQuadFast.h"
#include "operator/bilinearForm.hh"
#include "hp2D/bf_laplace.hh"

namespace concepts {
namespace gfem {

class GfemLaplace : 
  //public concepts::BilinearForm<Cmplx,Real> 
  public hp2D::Laplace<Cmplx> 
{
public:
    typedef RCP<const ElementFormula<Cmplx> >  ElementFormulaRCP;

    GfemLaplace(ElementFormulaRCP formula,
        bool ignoreNonGfem = false) 
      : hp2D::Laplace<Cmplx>(formula)
      , formula(formula)
      , ignoreNonGfem(ignoreNonGfem)
    { }


    virtual GfemLaplace* clone() const { 
      return new GfemLaplace(*this);
    }

    GfemLaplace(const GfemLaplace& other)
      : hp2D::Laplace<Cmplx>(other.formula)
      , formula(other.formula)
      , ignoreNonGfem(ignoreNonGfem)
    { }

    virtual void operator()(
          const concepts::Element<Real>& elmX,
          const concepts::Element<Real>& elmY, 
          concepts::ElementMatrix<Cmplx>& em); 

    void operator()(
        const GfemQuad& elmX, 
        const GfemQuad& elmY, 
        concepts::ElementMatrix<Cmplx>& em);

    void operator()(
        const GfemQuadFast& elmX, 
        const GfemQuadFast& elmY, 
        concepts::ElementMatrix<Cmplx>& em);

    void computeEM_reference(
          const GfemQuadFast& elmX, 
          const GfemQuadFast& elmY, 
          concepts::ElementMatrix<Cmplx>& em);


    virtual ~GfemLaplace() {}

protected:
    virtual std::ostream& info(std::ostream& os) const;

private:
    ElementFormulaRCP formula; 

protected:
    void allocHelpers(const GfemQuadFast& elm1);
    void fillHelpers(const GfemQuadFast& elm1);

    Cmplx& getH4(int ikappa1, int ikappa2, int imic1, int imic2) {
      assert(ikappa1 <  nKappa[0]);
      assert(ikappa2 <  nKappa[1]);
      assert(imic1 <  nMic);
      assert(imic2 <  nMic);
      return H4[ikappa1 + (ikappa2 + (imic1 + imic2 * nMic) * nKappa[1] ) * nKappa[0]];
    }

    Cmplx& getH2(int ikappa1, int imac2, int jmac2, int imic1, int imic2) {
      assert(ikappa1 <  nKappa[0]);
      assert(imac2 <  nMacro[1]);
      assert(jmac2 <  nMacro[1]);
      assert(imic1 <  nMic);
      assert(imic2 <  nMic);
      return H2[ikappa1 + (imac2 + (jmac2 + (imic1 + imic2 * nMic) 
                        * nMacro[1]) * nMacro[1] ) * nKappa[0]];
    }

    Cmplx& getH1(int nu1, int jmac1, int imac2, int jmac2, int imic1, int imic2) {
      assert(nu1 <  nMacro[0]);
      assert(jmac1 <  nMacro[0]);
      assert(imac2 <  nMacro[1]);
      assert(jmac2 <  nMacro[1]);
      assert(imic1 <  nMic);
      assert(imic2 <  nMic);
      return H1[nu1 + (jmac1 + (imac2 + (jmac2 + (imic1 + imic2 * nMic) 
                        * nMacro[1]) * nMacro[1] ) * nMacro[1]) * nMacro[1] ];
    }

protected:
    int nMic;
    int pMacro[2];
    int nMacro[2]; // = pMacro + 1
    int nKappa[2]; // = 2*pMacro + 1

    Cmplx* H4;
    Cmplx* H2;
    Cmplx* H1;

    static const int SELX = 0;
    static const int SELY = 1;
    static const int SELMAC = KarniadakisCoeffs::SELMAC;
    static const int SELMIC = 1 - KarniadakisCoeffs::SELMAC;
    
    int i_diff1; // differential selector func 1
    int i_diff2; // differential selector func 2
    int i_comp; // component selector

    KarniadakisCoeffs* kcoeffsX_;
    KarniadakisCoeffs* kcoeffsY_;
    bool ignoreNonGfem;
};

} // namespace
} // namespace