app-bholger/gfem/gfemQuad.h

Go to the documentation of this file.

#pragma once

#include <iostream>
#include <memory>
#include "basics/typedefs.hh"
#include "basics/outputOperator.hh"
#include "basics/vectorsMatrices.hh"
#include "geometry/cell2D.hh"
#include "geometry/integral.hh"
#include "integration/quadRule.hh"
#include "space/integral.hh"
#include "space/element.hh"
#include "function/vector.hh"
#include "integration/karniadakis.hh"
#include "hp2D/quad.hh"
#include "unitcell.h"
#include "gfemOverlapCondition.h"
#include "gfemDofIterator.h"
#include "geometry/cell2D.hh"

// If micro-oscillations are cutoff on the overlap edges, their corresponding edge and vertex
// dofs are disabled. If the micro-oscilations should still appear as "inner" dofs, this
// flag should be defined.
//#define USE_SPECIAL_OVERLAP 

//#define SPECIAL_MIC_ONE

//#define MIC_P_ONLY_FIRST

namespace concepts {
namespace gfem {


class GfemQuad : public concepts::ElementWithCell<Real> {
public:
  enum EvalMode { EVAL_ON_MICRO = 0, EVAL_ON_MACRO };

    GfemQuad(concepts::Quad2d& macroC, const Unitcell& uc
        , concepts::TColumn<Real>* T, int p_macro[2]
        , GfemOverlapCondition& overC, EvalMode eval_mode = EVAL_ON_MICRO); 
    virtual ~GfemQuad();

private:
    GfemQuad(const GfemQuad& other);
    GfemQuad& operator=(const GfemQuad& other);

public:

    virtual const Quad2d& cell() const {
      return macroC;
    }

    void recomputeShapefunctions() {
      clear();
      initialise();
    }
    
    static int disableFull(int mask, int i_edge) {
      return mask | (1024 << i_edge);
    }
 
    static int disablePartly(int mask, int i_elm) {
      return mask | (1 << i_elm);
    }
   
    static bool partlyDisabled(int mask, int i_elm) {
      return (mask & (1024 << i_elm)) | (mask & (1 << i_elm));
    }

    static bool fullyDisabled(int mask, int i_elm) {
      return mask & (1024 << i_elm);
    }
 
    static bool partlyDisabled(int mask) {
      return mask != 0;
    }
   
    static bool fullyDisabled(int mask) {
      return mask > 1024;
    }

    virtual const TMatrixBase<Real>& T() const {
      return T_; 
    }

    

    virtual const ElementGraphics<Real>* graphics() const;

    const concepts::Quad& support () const
    {
      // TODO: maybe use macroC.support()
      return macroC.connector();
    }

    MapReal2d getMacroJacobian(const Real& xi, const Real& eta) const {
      return macroC.jacobian(xi, eta);
    }

    MapReal2d getMacroJacobian(int i_quad) const {
      const Real pos_x = abscissas_x[i_quad];
      const Real pos_y = abscissas_x[i_quad];
      // scaling due to the fact that quadrature rules operate on [-1, 1]
      //return (getMacroJacobian((pos_x+1)/2, (pos_y+1)/2) *= .5);
      return (getMacroJacobian((pos_x+1)/2, (pos_y+1)/2));
    }

    int getNRepX() const {
      return n_repX;
    }

    int getNRepY() const {
      return n_repY;
    }

    MapReal2d getMicroJacobian(int i_quad) const {
      return (uc_jacobian[i_quad] * 2).scaleCols(Real2d(n_repX, n_repY));
    }

    const Cmplx& getMicroShape(int i_micro, int i_quad) const {
      return const_cast<GfemQuad*>(this)->uc_shape_val(i_micro, i_quad);
    }

    const Cmplx& uc_shape_val(int i_ucfun, int i_quad) const {
      return const_cast<GfemQuad*>(this)->uc_shape_val(i_ucfun, i_quad);
    }

    Cmplx& uc_shape_val(int i_ucfun, int i_quad) {
      assert(i_quad < nquad);
      assert(i_ucfun < numMicro());
      return uc_shape[i_ucfun * nquad + i_quad];
    }

    const Cmplx2d& uc_shape_valD(int i_ucfun, int i_quad) const {
      return const_cast<GfemQuad*>(this)->uc_shape_valD(i_ucfun, i_quad);
    }

    Cmplx2d& uc_shape_valD(int i_ucfun, int i_quad) {
      assert(i_quad < nquad);
      assert(i_ucfun < numMicro());
      return uc_shape_D[i_ucfun * nquad + i_quad];
    }

    const Real& getWeight(int i_quad) const  {
      return weights[i_quad];
    }

    Real2d getWorldPoint(int i_quad) const  {
      return macroC.chi((abscissas_x[i_quad]+1)/2, (abscissas_y[i_quad]+1)/2);
    }

    Real2d getUCPoint(int i_quad) const  {
      return Real2d((abscissas_x[i_quad]+1)/2, (abscissas_y[i_quad]+1)/2);
    }

    Real2d getAbscissa(int i_quad) const  {
      return Real2d(abscissas_x[i_quad], abscissas_y[i_quad]);
    }

    Real getMacroShapeX(int i_func, int i_quad) const  {
      if(i_func < 0)
        return 1;
      return macro_shapeX->values()[i_func*nquad + i_quad];
    }

    Real getMacroShapeY(int i_func, int i_quad) const  {
      if(i_func < 0)
        return 1;
      return macro_shapeY->values()[i_func*nquad + i_quad];
    }

    Real2d getMacroShapeD(int i_funcX, int i_funcY, int i_quad) const  {
      // factor 2 due to [-1,1] -> [0, 1]
      return Real2d(2*getMacroShapeDX(i_funcX, i_quad)*getMacroShapeY(i_funcY, i_quad)
          ,         2*getMacroShapeDY(i_funcY, i_quad)*getMacroShapeX(i_funcX, i_quad)  );
    }

    Real getMacroShapeDX(int i_func, int i_quad) const  {
      if(i_func < 0)
        return 0;
      return macro_shapeDX->values()[i_func*nquad + i_quad];
    }

    Real getMacroShapeDY(int i_func, int i_quad) const  {
      if(i_func < 0)
        return 0;
      return macro_shapeDY->values()[i_func*nquad + i_quad];
    }

    static int packLdofUC(int i_shape_x, int i_shape_y, int n_shape_macro_x) {
      assert(i_shape_x < n_shape_macro_x);
      return i_shape_y * n_shape_macro_x + i_shape_x;
    }

    int numMacro() const {
      return macro_shapeX->n() * macro_shapeY->n();
    }
    
    int numMacroX() const {
      return macro_shapeX->n();
    }

    int numMacroY() const {
      return macro_shapeY->n();
    }

    int numMicro() const {
      return n_micro;
    }

    int numUCFun() const {
      return n_ucfun;
    }

    int numQuadPoints() const {
      return nquad;
    }

    const Real* getAbscissas_x() const {
      return abscissas_x;
    }

    const Real* getAbscissas_y() const {
      return abscissas_y;
    }

    const Real* getWeights() const {
      return weights;
    }

    const Cmplx* getUCShape() const {
      return uc_shape;
    }

    const Cmplx2d* getUCShapeD() const {
      return uc_shape_D;
    }

    const Unitcell& getUC() const {
      return uc;
    }

    int num_shape() const {
      return numTidx();
    }

    const concepts::Karniadakis<1,0>& getMacroShapeX() const {
      return *macro_shapeX;
    }

    const concepts::Karniadakis<1,0>& getMacroShapeY() const {
      return *macro_shapeY;
    }

    const concepts::Karniadakis<1,1>& getMacroShapeDX() const {
      return *macro_shapeDX;
    }

    const concepts::Karniadakis<1,1>& getMacroShapeDY() const {
      return *macro_shapeDY;
    }

    const int* getPMax() {
      return p_max;
    }

    const int* getPMacro() {
      return p_macro;
    }

    const ElementWithCell<Real>* getMicroElem(int i_quad) const {
      return micro_elems[i_quad];
    }

    int numTidx() const {
      return numLDof();
    }

    int numLDof() const {
      return begin().numLdofs();
    }


    EvalMode getEvalMode() const {
      return eval_mode;
    }

    const GfemDofIterator& begin() const {
      return it_begin;
    }

protected:
    void initialise();
    void clear();

private:
    Quad2d& macroC;
    const Unitcell& uc;
    concepts::TMatrix<Real> T_;
    GfemOverlapCondition overC;
    GfemDofIterator it_begin;
    int p_max[2];
    int p_macro[2];
    int n_micro;
    int n_borderFuns;
    int n_ucfun; 

    int nquad;
    int n_shape_macro_x;
    int n_shape_macro_y;

    Real offX;
    Real offY;
    Real sizeX;
    Real sizeY;
    Real* abscissas_x;
    Real* abscissas_y;
    Real* weights;
    MapReal2d* uc_jacobian;

    Cmplx* uc_shape;
    Cmplx2d* uc_shape_D;
    const ElementWithCell<Real>** micro_elems;

    Karniadakis<1,0>* macro_shapeX;
    Karniadakis<1,0>* macro_shapeY;
    Karniadakis<1,1>* macro_shapeDX;
    Karniadakis<1,1>* macro_shapeDY;
  
    int n_repX;
    int n_repY;

    EvalMode eval_mode;
};

} // namespace
} // namespace