implicitEquilibratedAO.hh

Go to the documentation of this file.

 
#ifndef implicitAO_hh
#define implicitAO_hh
 
#include "fluxes.hh"
//#include "innerResidual.hh"
//#include "linInnerProd.hh"
//#include "moments.hh"
//#include "patches.hh"
 
#include "geometry/mutableMesh.hh"
//for computing L2Error
 
namespace concepts{
 
class ImplicitEquilibratedA0Error: public ExceptionBase {
public:
  ImplicitEquilibratedA0Error(const std::string& errorMessage) throw () {
    errorMessage_ = errorMessage;
  }
 
  virtual ~ImplicitEquilibratedA0Error() throw () {
  }
  ;
protected:
  virtual std::ostream& info(std::ostream& os) const throw () {
    return os << concepts::typeOf(*this) << "[" << errorMessage_ << " ]";
 
  }
private:
  std::string errorMessage_;
};
}
 
 
 
 
namespace hp2D{
 
//class representing the error estimator presented by Ainsworth and Oden in 2d
// based on equilibrated moments
 
template<class F>
class ImplicitEquilibratedA0: public concepts::LocalEstimator<F>{
 
public:
 
  ImplicitEquilibratedA0(const concepts::SpaceOnCells<F>& spc,
               const concepts::Vector<F>& sol,
               concepts::LinearForm<F>& f,
               ushort loc_rule = 2 ,
               const concepts::ElementFormula<F>* a = 0):
    concepts::LocalEstimator<F>(spc, sol), rule_(loc_rule), f_(f){
    //insert zero to inner edges attributes since by default inner edges have
    iAttrbs_.insert(0);
    concepts::Estimator<F>::globError_ = 0.0;
 
    //the implementation of  the function have to be done
    if(a)
      throw concepts::MissingFeature("Not implemented");
  };
 
 
  virtual ~ImplicitEquilibratedA0(){};
 
  void addLinearformInfo(concepts::LinearForm<F>& form, const F scalar =1){
    B_.push_back(&form);
    scalars_.push_back(scalar);
  }
 
  void addBillinearInfo(concepts::BilinearForm<F>& bform, const F scalar = 1){
    L_.push_back(&bform);
    scalars_L.push_back(scalar);
  }
 
  void addNeumannInfo(concepts::ParsedFormula<F>& frm, const concepts::Set<uint>& attrb){
    nfrms_.push_back(frm);
    nAttrbs_.push_back(attrb);
    for(concepts::Set<uint>::const_iterator iter = attrb.begin(); iter != attrb.end(); ++iter)
      help_[*iter] = nfrms_.size()-1;
  }
  void addDirichletInfo(const concepts::Set<uint>& attrb){
  //  concepts::Set<uint>::const_iterator iter = attrb.begin();
  //  for(;iter != attrb.end();++iter){
    //union dirichlet data
    dAttrbs_ = dAttrbs_||attrb;
  //  }
 
  }
 
  void addInnerInfo(const concepts::Set<uint>& attrb){
    iAttrbs_= iAttrbs_||attrb;
  }
 
 
  void isunique(bool unique=true){
    unique_=unique;
    //TODO: the action is setted
  }
 
  //Computes the error with all information, make sure to add information before with
  void compute(){
 
    const hp2D::hpAdaptiveSpaceH1* hpSpc =dynamic_cast<const hp2D::hpAdaptiveSpaceH1*> (&(this->spc_));
    if (hpSpc) {
 
      const concepts::BoundaryConditions* bcP = hpSpc->helper().bc();
      //copy to undo const
      concepts::BoundaryConditions bc(*bcP);
      //Patches
      geometry::VtxToPatchMaps patchMaps(*hpSpc, bc);
 
      //Build the inner Residual with the help of the added residual information
      //the Linearforms are just evaluated on boundary basis functions, therefore
      //the linearforms must be choiceable, i.e. a LinearFormChoice
      concepts::InnerResidual<F> innerRes(*hpSpc);
      concepts::LinearFormChoice* lfc=0;
      //add B(u_h,v)
      for(uint i = 0; i < B_.size(); ++i ){
        lfc = dynamic_cast<concepts::LinearFormChoice*>(B_[i]);
        if(lfc){
          lfc->setBasis(concepts::BND);
          innerRes.add(*(B_[i]), (-1)*scalars_[i]);
        }else{
         throw conceptsException(concepts::ImplicitEquilibratedA0Error("Linearform must be a LinearformChoice"));
        }
      }
      //add f(v)
      lfc = dynamic_cast<concepts::LinearFormChoice*>(&f_);
      if(lfc){
        lfc->setBasis(concepts::BND);
        innerRes.add(f_, -1);
      }else
         throw conceptsException(concepts::ImplicitEquilibratedA0Error("Linearform must be a LinearformChoice"));
 
 
 
      //Build approximated Moments
      hp2D::ApproxMoments<F> apprxMoments(*hpSpc, this->sol_);
      //add Neumann data
      for(uint i = 0 ; i < nfrms_.size(); ++i)
        apprxMoments.addNeumann(nfrms_[i],nAttrbs_[i] );
      //add Dirichlet data
      if(dAttrbs_.size())
        apprxMoments.addDirichlet(dAttrbs_);
 
      apprxMoments.addComplete();
 
      //compute equilibrated Moments
      hp2D::EquilibratedMomentsAO<F> moments(*hpSpc, patchMaps, innerRes, apprxMoments);
 
      //build the Fluxes for local neumannboundary condition
      hp2D::TraceSpace fluxSpc(*hpSpc,iAttrbs_);
      hp2D::Fluxes flux(fluxSpc, moments);
 
 
      //Reset Linearforms to evaluate on all basis functions, use for local problems
      for(uint i = 0; i < B_.size(); ++i ){
        //now lfc is ok, was tested before
        lfc = dynamic_cast<concepts::LinearFormChoice*>(B_[i]);
        lfc->setBasis(concepts::Default);
      }
      lfc = dynamic_cast<concepts::LinearFormChoice*>(&f_);
      lfc->setBasis(concepts::Default);
 
      //squared local errors
      concepts::HashMap<Real> error_2;
      error_2.rehash(hpSpc->nelm());
 
 
      //if unique local problem compute standard
      if(unique_) {
 
        std::unique_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(hpSpc->scan());
        while (!sc->eos()) {
          hp2D::Element<Real>& elm = (*sc)++;
          hp2D::Quad<Real>* quad = dynamic_cast<hp2D::Quad<Real>*> (&elm);
          if (quad) {
            //build the local Space
            concepts::MutableMesh2 quadMesh;
            concepts::Quad2d* quad2d =  dynamic_cast<concepts::Quad2d*>(&quad->cell());
            if (quad2d) {
              quadMesh.addCell(quad2d, false);
              uint p = std::max(quad->p()[0], quad->p()[1]);
              hp2D::hpAdaptiveSpaceH1 spc_K(quadMesh, 0, p+ rule_, &bc);
              spc_K.rebuild();
 
 
              //solution of the local problem
              //concepts::Vector<F> sol_K;
 
              //system Matrix
              concepts::SparseMatrix<F> A_K(spc_K.dim());
              //rhs b_K
              concepts::Vector<F> b_K(spc_K);
              //build the local right hand side, rhs = f(v) - B(u_h,v) + g(v) + flux(v)
              buildlocalRhs_(quad2d->connector(), spc_K, flux, bc, b_K);
 
              //build the systemmatrix
              buildlocalA_(spc_K, A_K);
 
              //solve the local problem
              concepts::Vector<F> sol_K(spc_K);
              //TODO: Have Mumps? if have mumps else superLU else gmres  TODO: use mumps for speed up
              concepts::SuperLU<F> solver_K(A_K);
              solver_K(b_K, sol_K);
 
 
              //compute the local error squared
              error_2[quad2d->connector().key()] =  computeLocalError_2(spc_K,sol_K);
 
            }//if quad 2d
            else
              throw concepts::MissingFeature("Not implemented A ");
          }//if quad
          else
            throw concepts::MissingFeature("Not implemented B ");
        }//while
 
      }// unique case
 
      std::unique_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(hpSpc->scan());
        while (!sc->eos()) {
          hp2D::Element<Real>& elm = (*sc)++;
          hp2D::Quad<Real>* quad =
              dynamic_cast<hp2D::Quad<Real>*> (&elm);
          if (quad) {
            //build the local Space
            concepts::MutableMesh2 quadMesh;
            concepts::Quad2d* quad2d =
                dynamic_cast<concepts::Quad2d*> (&quad->cell());
            if (quad2d) {
              quadMesh.addCell(quad2d, false);
              uint p = std::max(quad->p()[0], quad->p()[1]);
              hp2D::hpAdaptiveSpaceH1 spc_K(quadMesh, 0, p + rule_, &bc);
              spc_K.rebuild();
 
              //solution of the local problem
              //concepts::Vector<F> sol_K;
 
              //rhs b_K
              concepts::Vector<F> b_K(spc_K);
              //build the local right hand side, rhs = f(v) - B(u_h,v) + g(v) + flux(v)
              bool unique = buildlocalRhs_(quad2d->connector(), spc_K, flux, bc,b_K);
 
            //system Matrix
              concepts::SparseMatrix<F> A_K(spc_K.dim());
 
 
              //build the systemmatrix
              buildlocalA_(spc_K, A_K);
 
              //solve the local problem
              concepts::Vector<F> sol_K(spc_K);
              //TODO: Have Mumps? if have mumps else superLU else gmres  TODO: use mumps for speed up
              concepts::SuperLU<F> solver_K(A_K);
              solver_K(b_K, sol_K);
 
            }
          }//quad
        }//while
 
 
 
        //        hp2D::Riesz<Real> lf_one(concepts::ConstFormula<Real>(1.0));
        //        concepts::Vector<Real> v_one_K(spcOneElement, lf_one);
        //        uint dim_K = spcOneElement.dim();
        //        concepts::SparseMatrix<Real> A_K(dim_K + 1, dim_K + 1);
        //        concepts::IndexRange idxrange_K(0, dim_K - 1);
        //        concepts::Set<concepts::IndexRange> setIdx_K;
        //        setIdx_K.insert(idxrange_K);
        //        concepts::SubMatrixN<concepts::SparseMatrix<Real> > MsubA_K(
        //            A_K, setIdx_K, setIdx_K);
        //        MsubA_K.assembly(MsubA_K, spcOneElement, la_K);
        //        for (uint i = 0; i < dim_K; i++)
        //          A_K(i, dim_K) = A_K(dim_K, i) = v_one_K[i];
        //        A_K.compress();
        //
        //
        //
        //        concepts::Vector<Real> bigrhs_K(dim_K+1);
        //        for(uint i = 0 ; i < dim_K; i++)
        //          bigrhs_K[i]=rhs_K[i];
        //        bigrhs_K[dim_K]=0;
        //
        //
        //        //solve the System
        //        concepts::SuperLU<Real> solver_K(A_K);
        //        //concepts::CG<Real> solver(A, 1e-6, 1000);
        //        concepts::Vector<Real> bigsol_K(dim_K + 1);
        //        solver_K(bigrhs_K, bigsol_K);
        //        sol_K.resize(dim_K);
        //        for(uint i = 0 ; i < dim_K ; i++)
        //          sol_K[i]= bigsol_K[i];
        //
        //        //sol_K(bigsol_K, setIdx_K);
 
 
 
 
 
 
 
      concepts::LocalEstimator<F>::locError_.rehash(hpSpc->nelm());
 
      //Place local errors and global error
 
 
 
      for(concepts::HashMap<Real>::const_iterator iter_err = error_2.begin(); iter_err != error_2.end(); ++iter_err){
        concepts::LocalEstimator<F>::locError_[iter_err->first] = sqrt(iter_err->second);
        concepts::Estimator<F>::globError_ += iter_err->second;
      }
      concepts::Estimator<F>::globError_ = sqrt(concepts::Estimator<F>::globError_);
 
    }//if hpSpc
    else
      throw conceptsException(concepts::MissingFeature("Just AdaptiveSpace implemented."));
 
 
 
 
 
  }
 
  virtual ImplicitEquilibratedA0* clone() const{return new ImplicitEquilibratedA0(*this);}
 
protected:
 
  virtual std::ostream& info(std::ostream& os) const{
    return os<< concepts::typeOf(*this);
  }
 
private:
 
  //weak residual as informational part
  concepts::Sequence<concepts::LinearForm<F>* > B_;
  //scalar multiplication factors for the wres_ linearforms
  concepts::Sequence<F> scalars_;
  //rhs lform
  concepts::LinearForm<F>& f_;
 
  //Billinearinformation
  concepts::Sequence<concepts::BilinearForm<F>* > L_;
  concepts::Sequence<F> scalars_L;
 
 
 
  //flag mark if local problems are solveable uniquely,
  //i.e. for pure laplace problems local problems are solveable but not unique, since a local
  //problem is a pure NeumannProblem
  bool unique_;
  //local polynomial addition rule, elementpolynomial.max + rule_
  ushort rule_;
 
  //Neumann formulas
  concepts::Sequence<concepts::ParsedFormula<F> > nfrms_;
  //Neumann Attributes
  concepts::Sequence<concepts::Set<uint> > nAttrbs_;
  //help structur to find position of neumann formula
  concepts::HashMap<uint> help_;
  //Dirichlet Attributes
  concepts::Set<uint> dAttrbs_;
  //inner Edges Attributes
  concepts::Set<uint> iAttrbs_;
 
 
 
  //TODO: Make this const
  //build rhs for the unique case
  bool buildlocalRhs_(const concepts::Quad& quad,
            const hp2D::hpAdaptiveSpaceH1& spc_K,
            const hp2D::Fluxes& flux,
            const concepts::BoundaryConditions& bc,
            concepts::Vector<F>& rhs){
 
 
    // build rhs to int flux(v) ds on interior edges
    hp2D::TraceSpace tspcFluxElements(spc_K, iAttrbs_);
 
    //get all possible neumann edges
    concepts::Set<uint> nattrb_K;
    for (uint i = 0; i < 4; ++i) {
      if (bc(quad.edge(i)->attrib()).type()
          == concepts::Boundary::NEUMANN)
        nattrb_K.insert(quad.edge(i)->attrib().attrib());
    }
 
 
    //number of pure Neumannedges, it is 4 if it's an interior quad or a quad that intersects neumannbdn
    uint noN = tspcFluxElements.nelm() + nattrb_K.size();
    //up to 4 edges can be neumann ones
    conceptsAssert(noN<5, concepts::Assertion());
 
    //is a pure non unique neumannproblem
 
    //local problem is unique solveable or touches a boundary element that is not Neumann (i.e. a Dirichlet)
    //TODO: What if different boundary type?
    if (unique_ || noN < 4){
 
      hp2D::LocalFluxes locFlux(tspcFluxElements, quad.key(), flux);
      hp1D::Riesz<F> lformFluxes_new(locFlux);
      //build rhs vector coressponding to the fluxes
      concepts::Vector<F> rhs_flux(tspcFluxElements, lformFluxes_new);
 
      // f_K(v)
      concepts::Vector<F> rhs_f_K(spc_K, f_);
 
      //initialisation plus data copy
      rhs = rhs_flux;
      rhs += rhs_f_K;
 
      //TODO: was ist wenn die billinearform nicht auf dem Space lebt?
      for (uint i = 0; i < B_.size(); ++i) {
        concepts::Vector<F> rhs_B(spc_K, *(B_[i]));
        rhs += rhs_B * scalars_[i];
      }
 
      //add neumann contributions
      if (nattrb_K.size() > 0) {
        for (concepts::Set<uint>::const_iterator iter =
            nattrb_K.begin(); iter != nattrb_K.end(); ++iter) {
          hp2D::TraceSpace tspcNeumann(spc_K, *iter);
          hp1D::Riesz<F> lform_g(nfrms_[help_[*iter]]);
          //Boundary g_K integrals
          concepts::Vector<F> rhs_g(tspcNeumann, lform_g);
          rhs += rhs_g;
        }
      }
 
      //problem is unique solveable return true
      return true;
 
    }else{
      //local problem has solution that is unique up to a constant
 
 
 
 
 
    }
 
 
 
 
 
 
 
 
  }//build local Rhs
 
 
 
 
 
  //TODO: Make this const
  //build rhs for the unique case
  void buildnonUniquelocalRhs_(const concepts::Quad& quad,
            const hp2D::hpAdaptiveSpaceH1& spc_K,
            const hp2D::Fluxes& flux,
            const concepts::BoundaryConditions& bc,
            concepts::Vector<F>& rhs){
 
 
    // build rhs to int flux(v) ds on interior edges
    hp2D::TraceSpace tspcFluxElements(spc_K, iAttrbs_);
 
    //get all possible neumann edges
    concepts::Set<uint> nattrb_K;
    for(uint i=0; i< 4 ; ++i){
      if(bc(quad.edge(i)->attrib()).type()== concepts::Boundary::NEUMANN)
        nattrb_K.insert(quad.edge(i)->attrib().attrib());
    }
 
 
 
    hp2D::LocalFluxes locFlux(tspcFluxElements, quad.key(), flux);
    hp1D::Riesz<F> lformFluxes_new(locFlux);
    //build rhs vector coressponding to the fluxes
    concepts::Vector<F> rhs_flux(tspcFluxElements, lformFluxes_new);
 
    // f_K(v)
    concepts::Vector<F> rhs_f_K(spc_K, f_);
 
    //initialisation plus data copy
    rhs = rhs_flux;
    rhs+= rhs_f_K;
 
    //TODO: was ist wenn die billinearform nicht auf dem Space lebt?
    for(uint i = 0 ; i < B_.size(); ++i){
      concepts::Vector<F> rhs_B(spc_K,*(B_[i]));
      rhs+= rhs_B*scalars_[i];
    }
 
 
    //add neumann contributions
    if(nattrb_K.size()>0){
      for(concepts::Set<uint>::const_iterator iter = nattrb_K.begin(); iter != nattrb_K.end(); ++iter){
        hp2D::TraceSpace tspcNeumann(spc_K, *iter);
        hp1D::Riesz<F> lform_g(nfrms_[help_[*iter]]);
        //Boundary g_K integrals
        concepts::Vector<F> rhs_g(tspcNeumann, lform_g);
        rhs += rhs_g;
      }
    }
  }//build local Rhs
 
 
 
 
 
 
 
 
 
 
 
 
 
  void buildlocalA_(const hp2D::hpAdaptiveSpaceH1& spc_K, concepts::SparseMatrix<F>& A_K){
    typename concepts::Sequence<concepts::BilinearForm<F>* >::const_iterator iter =  L_.begin();
    typename concepts::Sequence<F>::const_iterator iter_s =  scalars_L.begin();
    //Billinearinformation
    for(; iter!=L_.end(); ++iter, ++iter_s){
      concepts::SparseMatrix<F> A(spc_K, **iter);
      //add local contribution to the whole local matrix
      A.addInto(A_K,*iter_s);
    }
  }
 
 
  Real computeLocalError_2(hp2D::hpAdaptiveSpaceH1& spc_K, const concepts::Vector<F>& sol_K){
    //TODO: Write this general !!
    concepts::ElementFormulaVector<1,F> errID_K(spc_K, sol_K, hp2D::Value<Real>());
    concepts::ElementFormulaVector<2,F> errGrad_K(spc_K, sol_K, hp2D::Grad<Real>());
    if(unique_){
    Real grad_l2 = concepts::L2product(spc_K,errGrad_K);
    Real id_l2 = concepts::L2product(spc_K, errID_K);
    return  grad_l2 + id_l2;
    }else{
      //H1 seminorm
      return concepts::L2product(spc_K,errGrad_K);
    }
  }
 
 
 
 
 
};
 
 
} //namespace
 
#endif // implicitAO_hh