app-bholger/gfem/gfemIdentity.h

Go to the documentation of this file.

#pragma once

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

namespace concepts {
namespace gfem {

class GfemIdentity : public hp2D::Identity<Cmplx>
{
public:
    typedef concepts::RCP<const concepts::ElementFormula<Cmplx> >
      ElementFormulaRCP;

    GfemIdentity(ElementFormulaRCP formula) 
      : hp2D::Identity<Cmplx>(formula)
    { }

    GfemIdentity(const GfemIdentity& other)
      : hp2D::Identity<Cmplx>(other)
    { }


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

    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 ~GfemIdentity() {}

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

protected:
    void allocHelpers(const GfemDofIteratorFast& itBeg);
    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;
};

} // namespace
} // namespace