app-bholger/gfem/gfemSpace.h

Go to the documentation of this file.

#pragma once

#include "gfemQuad.h"
#include "space/space.hh"
#include "geometry/mesh.hh"
#include "hp2D/quad.hh"
#include <vector>
#include <map>
#include <iostream>
#include "limits"
#include "gfemControl.h"

namespace concepts {
namespace gfem {

struct GfemSpaceParams {
  GfemSpaceParams(Attribute attr_microOsc,  int p_macro = 1, int p_outer = 1, int p_edge = 1e6)
  {
    addGfemAttr(attr_microOsc, p_macro);

    p_outers[0] = p_outer;
    p_outers[1] = p_outer;

    p_edges[0] = std::min(p_outers[0], p_edge); 
    p_edges[1] = std::min(p_outers[1], p_edge); 
 }

  GfemSpaceParams(int p_outer = 1, int p_edge = 1e6)
  {
    p_outers[0] = p_outer;
    p_outers[1] = p_outer;
    p_edges[0] = std::min(p_outers[0], p_edge); 
    p_edges[1] = std::min(p_outers[1], p_edge); 
  }

  GfemSpaceParams(int po[2], int p_edge = 1e6)
  {
    p_outers[0] = po[0];
    p_outers[1] = po[1];
    p_edges[0] = std::min(p_outers[0], p_edge); 
    p_edges[1] = std::min(p_outers[1], p_edge); 
  }

  void addGfemAttr(Attribute gfemAttr, int p_macro) {
#if 0
    printf("****** FIXME: hardcoded barrier geometry! ****\n");
    if(int(gfemAttr) == 8) {
      gfemParamMap[gfemAttr] = GfemCellParams(p_macro
          , GfemCellParams::SEL_ALL, GfemCellParams::SEL_POS);
    } else if(int(gfemAttr) == 9) {
      gfemParamMap[gfemAttr] = GfemCellParams(p_macro
          , GfemCellParams::SEL_ALL, GfemCellParams::SEL_PER);
    } else {
      gfemParamMap[gfemAttr] = GfemCellParams(p_macro);
    }
#else
    gfemParamMap[gfemAttr] = GfemCellParams(p_macro);
#endif
  }

  const GfemCellParams* getCellParams(Attribute gfemAttr)  const
  {
    GfemParamMap::const_iterator it = gfemParamMap.find(gfemAttr);
    if(it == gfemParamMap.end())
      return NULL;
    return &(it->second);
  }

  int p_outers[2];
  int p_edges[2];

  typedef std::map<Attribute, GfemCellParams> GfemParamMap;
  GfemParamMap gfemParamMap; 
};


class GfemSpace : public concepts::SpaceOnCells<Real> 
{
public:
  typedef hp2D::Quad<Real> hpQuad;
  typedef StlVectorScan<ElementWithCell<Real> > Scan;
  typedef StlVectorScan<GfemQuad > ScanGfem;
  typedef StlVectorScan<hpQuad > ScanOrdinary;

  const static bool NO_SCALE_RES = GfemControl<0>::NO_SCALE_RES;
  const static bool USE_SCALE_RES = GfemControl<0>::USE_SCALE_RES;


  GfemSpace(Mesh2& msh, const Unitcell& uc, const GfemSpaceParams& params,
      const concepts::BoundaryConditions& bc);

  virtual ~GfemSpace();

  void rebuild();

  virtual uint dim() const {
    return dim_; 
  }

  int getDim() {
    return dim();
  }

  int getDimGfem() const {
    return dimGfem;
  }

  int getDimOverl() const {
    return dimOverl;
  }

  int genVertDof(GfemControlVert ctrl);

  virtual uint nelm() const {
    return elems.size();    
  }

  void clear();

  int getNelm() const {
    return elems.size();
  }

  int getNelmGfem() const {
    return elemGfemSubs.size(); 
  }

  int getNelmOrdinary() const {
    return elemGfemSubs.size();  
  }

  virtual Scan* scan() const {
    return new Scan(elems, elems.begin());
  }


  ScanOrdinary* scanOrdinary() const {
    return new ScanOrdinary(elemOrdinarySubs, elemOrdinarySubs.begin());
  }

  ScanGfem* scanGfem() const {
    return new ScanGfem(elemGfemSubs, elemGfemSubs.begin());
  }


  void assignVertexDofs(const concepts::Quad& cell_t, const GfemOverlapCondition& overc,
        concepts::TColumn<Real>*& T1);

  void assignEdgeDofs(const concepts::Quad& cell_t, const GfemOverlapCondition& overc,
        concepts::TColumn<Real>*& T1);

  void assignInnerDofs(const concepts::Quad& cell_t, const GfemOverlapCondition& overc,
        concepts::TColumn<Real>*& T1);

  bool isVertexDof(const concepts::Connector2& cntr, uint i) const;

  bool isEdgeDof(const concepts::Connector2& cntr, uint i) const;

  bool isInnerDof(const concepts::Connector2& quad) const;

  void verbosePrint();

  void recomputeShapefunctions();

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

  int allocDofs(int numDofs);

private:

  void processMesh(const Quad2d& quad);
  
  void processMicro(const Quad2d& quad);

  void processOverlap(const Quad2d& cell);

  void processBC(const Quad2d& quad);

  concepts::Mesh2& mesh;

  std::vector< ElementWithCell<Real>* > elems; 
  std::vector< GfemQuad* > elemGfemSubs; 
  std::vector< hpQuad* > elemOrdinarySubs; 
  typedef std::map< int, GfemControl<0> > ControlVert;
  typedef std::map< int, GfemControl<1> > ControlEdge;
  typedef std::map< int, GfemControl<2> > ControlCell;
  ControlVert controlVert;  
  ControlEdge controlEdge;  
  ControlCell controlCell;  

  int n_elem;
  int n_gfem;
  int n_ordinary;
  int dim_;
  int dimGfem;
  int dimOverl;

  GfemSpaceParams params;
  BoundaryConditions bc;
  const Unitcell& uc;

  //TODO: work here map of control informations
};

} // namespace
} // namespace