app-bholger/gfem/unitcell.h

Go to the documentation of this file.

#pragma once
#include "hp2D/space.hh"
#include "space/space.hh"
#include "vector"
#include "stdio.h"
#include "function/vector.hh"
#include "hp2D/hpAdaptiveSpaceH1.hh"
#include "toolbox/sharedPointer.hh"

namespace concepts {
namespace gfem {
  using concepts::Real;
  using concepts::Cmplx;

struct UC_Sol {
  typedef concepts::Vector<Cmplx> SolVec;

  UC_Sol(const concepts::Space<Real>& ucSpace, Cmplx k1=Cmplx(0,0), Cmplx k2=Cmplx(0,0))
    : vec(ucSpace)
    , sigma(0)
  {
    k[0] = k1;
    k[1] = k2;
  }

  UC_Sol(const SolVec& vec, Cmplx k1, Cmplx k2)
    : vec(vec)
    , sigma(0)
  { 
    k[0] = k1;
    k[1] = k2;
  }

  bool operator<(const UC_Sol& other) const {
    return sigma < other.sigma;
  }

  void print() const {
    printf("UC_Sol(%7.2f  %+7.2fi, %7.2f  %+7.2fi)\n"
        , real(k[0])
        , imag(k[0])
        , real(k[1])
        , imag(k[1])
        );
  }
  

  SolVec vec; // periodic coefficients of solution
  Cmplx k[2]; // k vector for quasi-periodic boundary conditions
  double sigma; // singular value, used in GfemBasisOne
};

struct GfemCellParams {
  static const int SEL_PER = -100; // only special mic osc (1) selected
  static const int SEL_ALL = 0;
  static const int SEL_NEG = -1;
  static const int SEL_POS = 1;

  GfemCellParams(int p_macro=0, int selX=SEL_ALL, int selY = SEL_ALL) {
    p_macros[0] = p_macro;
    p_macros[1] = p_macro;

    #if 0
    dir_sel[0] = SEL_ALL;
    dir_sel[1] = SEL_POS; // FIXME: this hardcodes infBarrier geometry
    //dir_sel[1] = SEL_ALL;
    #endif
    
    dir_sel[0] = selX;
    dir_sel[1] = selY;

  }

  int p_macros[2];

  int dir_sel[2];
};


struct Unitcell {
  Unitcell(RCP<const hp2D::hpAdaptiveSpaceH1> ucSpace);

  void addSolution(UC_Sol* sol) {
    uc_sols.push_back(sol);
  }

  const UC_Sol& getSolution(int i) const {
    return *(uc_sols[i]);
  }
  
  UC_Sol& getSolution(int i) {
    return *(uc_sols[i]);
  }

  int size() const {
    return (int)uc_sols.size();
  }


  ~Unitcell() {
    clear();
  }

  void clear() {
    if(!shallow) {
      for(int i=0; i < (int)uc_sols.size(); ++i)
        delete uc_sols[i];
    }
    uc_sols.clear();
  }

  Unitcell(const Unitcell& other, bool shallow = false)
    : ucSpace(other.ucSpace)
    , offX(other.offX)
    , offY(other.offY)
    , sizeX(other.sizeX)
    , sizeY(other.sizeY)
    , shallow(shallow)
  {
    if(!shallow) {
      // deep copy
      for(int i=0; i < (int)other.uc_sols.size(); ++i)
        uc_sols.push_back(new UC_Sol(*other.uc_sols[i]));
    } else {
      // shallow copy
      for(int i=0; i < (int)other.uc_sols.size(); ++i)
        uc_sols.push_back( other.uc_sols[i] );
    }
  }
private:
  Unitcell& operator=(const Unitcell& other);

public:

  RCP<const hp2D::hpAdaptiveSpaceH1> ucSpace;
  Real offX;
  Real offY;
  Real sizeX;
  Real sizeY;

private:
  std::vector< UC_Sol* > uc_sols;
  bool shallow;
};

struct UnitcellMask : public Unitcell {
  static const int UC_ACTIVE = 1;
  static const int UC_INACTIVE = 0;

  UnitcellMask(const Unitcell& other)
   : Unitcell(other, true)
   , uc_active(other.size(), UC_ACTIVE)
  {
    assert(size() == other.size());
    assert(sizeActive() == other.size());
    //new(*this) (other, false);
  }

  ~UnitcellMask() {
    uc_active.clear();
  }

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

public:
  static bool isActive(Cmplx k, int sel_mode) {
    if(sel_mode == GfemCellParams::SEL_ALL)
      return true;

    if(sel_mode == GfemCellParams::SEL_POS)
      return imag(k) >= 0;
    
    if(sel_mode == GfemCellParams::SEL_NEG)
      return imag(k) <= 0;

    if(sel_mode == GfemCellParams::SEL_PER)
      return imag(k) == 0 && real(k) == 0;

    if(sel_mode == GfemCellParams::SEL_ALL)
      return true;

    assert(false);
  }

  void updateActive(int i_sol, const GfemCellParams* cell_p) {
    setActive(i_sol, isActive(getSolution(i_sol), cell_p) );
  }

  void checkDeactivateAll(const GfemCellParams* cell_p) {
    for(int i=0; i < size(); ++i)
      checkDeactivate(i, cell_p);
  }

  void checkDeactivate(int i_sol, const GfemCellParams* cell_p) {
    if(!isActive(getSolution(i_sol), cell_p))
      setActive(i_sol, false );
  }

  static bool isActive(const UC_Sol& sol, const GfemCellParams* cell_p) {
    return isActive(sol.k[0], cell_p->dir_sel[0])
          && isActive(sol.k[1], cell_p->dir_sel[1]);
  }

#if 0
  const UC_Sol& getSolution(int i_sel, const GfemCellParams* cell_p) const {
    if(cell_p == NULL)
      return getSolution(i_sel);

    for(int i=0; i < size(); ++i) {
      const UC_Sol& sol = getSolution(i);

      if(isActive(sol, cell_p))
      {
        if(i_sel == 0)
          return sol;
        --i_sel;
      } 
    }
    assert(false);
  }


  int size(const GfemCellParams* cell_p) const {
    if(cell_p == NULL)
      return size();

    int active = 0;
    for(int i=0; i < size(); ++i) {
      const UC_Sol& sol = getSolution(i);

      if(isActive(sol, cell_p)) {
        active++;
      } else {
        printf("disabling ");
        sol.print();
      } 
    }
    return active;
  }
#endif
  bool isActive(int i_sol) const {
    return uc_active[i_sol] == UC_ACTIVE;
  }

  void setActive(int i_sol, bool state) {
    uc_active[i_sol] = state;
  }

  int sizeActive() const {
    assert(int(uc_active.size()) == size());

    int res = 0;
    for(int i=0; i < size(); ++i) {
      if(isActive(i))
        res++;
    }
    return res;
  }

private:
  std::vector< int > uc_active;
};

} // namespace
} // namespace