app-bholger/gfem/gfemUCQuad.h

Go to the documentation of this file.

#pragma once

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

namespace concepts {
namespace gfem {


class GfemUCQuad {
public:
  GfemUCQuad(const Unitcell& uc, int p_macro[2]
      ); 
  virtual ~GfemUCQuad();

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

public:
  void recomputeShapefunctions() {
    clear();
    initialise();

    std::vector< RCP<const ElementFormula<Cmplx> > > formulasOld(formulas);
    std::vector< EvalMode > modesOld(modes);
    formulas.clear();
    modes.clear();

    for(uint i=0; i < formulas.size(); ++i) {
      addFormula(formulasOld[i], modes[i]);
    }
  }
  
  int getFomulaId(const ElementFormula<Cmplx>* formula_adr) 
  {
    for(size_t i=0; i < formulas.size(); ++i) {
      if(formulas[i].get() == formula_adr)
        return i;
    }
    return -1;
  }

  void unsetCurFomula() {
    ucBF_lap = NULL;
    ucBF_id = NULL;
  }

  void setCurFomula(int id) {
    assert(id < (int)formulas.size());

    ucBF_id = ucBF_ids[id];
    ucBF_lap = ucBF_laps[id];
    cur_mode = modes[id];
    /*for(int i=0; i < N_LAP_COMP; ++i) { 
      ucBF_lap[i] = ucBF_laps[id][i];
    }*/
  }
  
  void activateFormula(const concepts::ElementFormula<Cmplx>* formula) {
    int id = getFomulaId(formula);
    if(id < 0) {
      std::cout << *formula << std::endl; 
      CONCEPTS_SIMPLE_EXC("formula not initialised ");
    }
    setCurFomula(id);
  }

  enum EvalMode { EVAL_ON_MICRO = 0, EVAL_ON_MACRO };

  void registerFormula(const RCP<const concepts::ElementFormula<Cmplx> >& formula, 
                       EvalMode mode = EVAL_ON_MICRO) 
  {
    int id = getFomulaId(formula.get());
    if(id == -1) {
      addFormula(formula, mode);

      //assert(getFomulaId(formula.get()) >= 0);
    } else {
      setCurFomula(id);
    }
  }

  int numMacro() const {
    return numMacroX() * numMacroY();
  }

  int numMacroX() const {
    return n_shape_macro_x;
  }

  int numMacroY() const {
    return n_shape_macro_y;
  }

  int numMicro() const {
    return n_micro;
  }

  int numUCFun() const {
    return n_ucfun;
  }

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

  const Cmplx* getQPpart(int iMic) const {
    return uc.getSolution(iMic).k;
  }

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

  const int* getPMax() {
    return p_max;
  }

  const int* getPMacro() {
    return p_macro;
  }

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

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

public: // access quadrature helper functions
  int numQuadPoints() const {
    return nquad;
  }

  Real2d getUCPoint(int i_quad, int repx=0, int repy=0) const  {
    assert(i_quad < numQuadPoints());
    return Real2d( ((abscissas_x[i_quad]+1)/2 + uc.offX + repx) * uc.sizeX 
        ,          ((abscissas_y[i_quad]+1)/2 + uc.offY + repy) * uc.sizeY );
  }

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

  MapReal2d getMicroJacobian(int i_quad) const {
    // TODO: check if we want it this way
    return (uc_jacobian[i_quad] * 2);
  }

    const Cmplx& getMicroShape(int i_micro, int i_quad) const {
      return const_cast<GfemUCQuad*>(this)->getMicroShape(i_micro, i_quad);
    }
    
    Real getMacroShape(int i_macro, Real x) const 
    {
      KarniadakisCoeffs& kcoeffs = KarniadakisCoeffs::getInstance();
      const PolyCoeff& p = kcoeffs.getCoeffs(i_macro); 
      return p.eval(x);
    }

    Real2d getMacroShapeD_unit(int i_macro1, int i_macro2, Real x, Real y) const 
    {
      KarniadakisCoeffs& kcoeffs = KarniadakisCoeffs::getInstance();
      const PolyCoeff& px   = kcoeffs.getCoeffs (i_macro1); 
      const PolyCoeff& py   = kcoeffs.getCoeffs (i_macro2); 
      const PolyCoeff& px_D = kcoeffs.getCoeffsD(i_macro1); 
      const PolyCoeff& py_D = kcoeffs.getCoeffsD(i_macro2); 
 
      const Real res_x = px_D.eval(x) * py  .eval(y);
      const Real res_y = px  .eval(x) * py_D.eval(y);

      return Real2d(res_x, res_y);
    }

    Real powi(Real x, int power) {
      Real res(1);
      for(int i = 0; i < power; ++i) {
        res *= x;
      }
      return res;
    }

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

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

    Cmplx2d& getMicroShapeD(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];
    }

    int numMonomX() {
      return 2*getPMax()[0]+1;
    }

    int numMonomY() {
      return 2*getPMax()[1]+1;
    }

    Cmplx getUCBF_Id(int px, int py, int iMicro1, int iMicro2) const {
      return const_cast<GfemUCQuad*>(this)->getUCBF_Id(px, py, iMicro1, iMicro2);
    }

    Cmplx& getUCBF_Id(int px, int py, int iMicro1, int iMicro2) {
      return getUCBF(ucBF_id, px, py, iMicro1, iMicro2);
    }
    Cmplx getUCBF_Lap(int px, int py, int iMicro1, int iMicro2,
        int i_diff1, int i_diff2, int i_comp) const 
    { 
      return const_cast<GfemUCQuad*>(this)->getUCBF_Lap(px, py, iMicro1, iMicro2,
          i_diff1, i_diff2, i_comp);
    }

    Cmplx& getUCBF_Lap(int px, int py, int iMicro1, int iMicro2, 
        int i_diff1, int i_diff2, int i_comp) 
    {
      Cmplx* arr = getUCBF_LapType(i_diff1, i_diff2, i_comp);
      return getUCBF(arr, px, py, iMicro1, iMicro2);
    }
protected:
    void initialise();
    void clear();
    void fillUCBF(const concepts::ElementFormula<Cmplx>& formula);

    Cmplx& getUCBF(Cmplx* vals,  int px, int py, int iMicro1, int iMicro2) {
      const int mul1 = numMicro();
      const int mul2 = mul1 * numMicro();
      const int mul3 = mul2 * numMonomX();

      return vals[iMicro1 + iMicro2*mul1 + px*mul2 + py*mul3];
    }

    Cmplx*& getUCBF_LapType(int i_diff1, int i_diff2, int i_comp) {
      const int idxType = i_diff1*4 + i_diff2*2 + i_comp;

      return ucBF_lap[idxType];
    }

    void addFormula(const RCP<const concepts::ElementFormula<Cmplx> >& formula,
                    EvalMode mode) 
    {
      formulas.push_back(formula);
      modes.push_back(mode);
      ucBF_ids.push_back(new Cmplx[n_ucBF]);
      ucBF_laps.push_back(new Cmplx*[8]);

      setCurFomula(formulas.size()-1);
      //ucBF_id = ucBF_ids.back();
      //ucBF_laps = ucBF_laps.back();

      getUCBF_LapType(0, 0, 0) = ucBF_id; // all derivatives are on macro functions 
      getUCBF_LapType(0, 0, 1) = ucBF_id; // all derivatives are on macro functions
      getUCBF_LapType(0, 1, 0) = new Cmplx[n_ucBF];
      getUCBF_LapType(0, 1, 1) = new Cmplx[n_ucBF];
      getUCBF_LapType(1, 0, 0) = new Cmplx[n_ucBF];
      getUCBF_LapType(1, 0, 1) = new Cmplx[n_ucBF];
      getUCBF_LapType(1, 1, 0) = new Cmplx[n_ucBF];
      getUCBF_LapType(1, 1, 1) = new Cmplx[n_ucBF];

      fillUCBF(*formula);
    }


private:
    const Unitcell& uc;

    GfemOverlapCondition overC;
    GfemDofIterator it_begin;

    std::vector< RCP<const ElementFormula<Cmplx> > > formulas;
    std::vector< EvalMode > modes;

    int p_max[2];
    int p_macro[2];
    int n_micro;
    int n_ucfun; 
    int n_ucBF;
    EvalMode cur_mode;

    int nquad;
    int n_shape_macro_x;
    int n_shape_macro_y;

    Real* abscissas_x;
    Real* abscissas_y;
    Real* weights;
    MapReal2d* uc_jacobian;
    const ElementWithCell<Real>** micro_elems;

    Cmplx* uc_shape;
    Cmplx2d* uc_shape_D;

    static const int N_LAP_COMP = 8;
    Cmplx* ucBF_id;
    Cmplx** ucBF_lap; // Cmplx*[N_LAP_COMP]

    std::vector<Cmplx*> ucBF_ids;
    std::vector<Cmplx**> ucBF_laps; // Cmplx*[N_LAP_COMP]
};

} // namespace
} // namespace