implicitResidual.hh

#ifndef implicitResidual_hh
#define implicitResidual_hh
 
#include "estimator.hh"
 
//TODO: Precise the includes
#include "hp2D.hh"
#include "formula.hh"
#include "basics/exceptions.hh"
 
#define COMPUTE_D 0
 
using concepts::Real;
 
namespace estimator{
 
 
 
template<class F>
class ImplicitResidual: public LocalEstimator<F>{
 
public:
 
 
 
  virtual ImplicitResidual<F>* clone() const = 0;
 
  void addBoundaryData(const concepts::Set<uint>& nSet,
             const concepts::ElementFormula<F>& g)
  {
  //  DEBUGL(1, " add boundary data");
      g_.push_back(&g);
      gAttrbs_.push_back(nSet);
 
      //eAttrbs_  =   eAttrbs_ || nSet;
      //DEBUGL(1, "all used attributes : " << eAttrbs_);
    //DEBUGL(1, " done... ");
  }
 
 
 
 
 
  void addBoundaryData(const concepts::Set<uint>& rSet,
             const concepts::ElementFormula<F>& h1,
             const concepts::ElementFormula<F>& h2)
  {
    throw conceptsException(concepts::MissingFeature("Not Supported"));
  }
 
 
 
protected:
 
 
  ImplicitResidual(const concepts::SpaceOnCells<F>& spc,
           const concepts::Vector<F>& sol):
      LocalEstimator<F>(spc, sol){};
 
  virtual std::ostream& info(std::ostream& os) const{
    return os << "Implicit Residual Error Estimator()";
  }
 
  //structure for the billinearforms of the local problems, with multiplikatives constant
  //concepts::Sequence< concepts::BilinearForm<F>* > bforms_;
  //const coeffs of the bforms
  //concepts::Sequence< F > c_bforms_;
 
  //local reference to the residual
  //const concepts::ElementFormula<F>& res_;
 
  //Neumann formulas
  concepts::Sequence<const concepts::ElementFormula<F>* > g_;
  //Neumann Attributes
  concepts::Sequence<concepts::Set<uint> > gAttrbs_;
 
 
  //all used edge attributes on that error estimates are placed.
  concepts::Set<uint> eAttrbs_;
 
  //map from edge attribute to the corresponding boundary formula, i.e. the neumann rhs
  //concepts::HashMap<const concepts::ElementFormula<F>* > g_;
 
 
};
 
 
 
 
 
 
 
 
 
 
 
 
template<class F>
class ImplicitResidual2D: public ImplicitResidual<F> {
 
 
 
public:
 
 
  ImplicitResidual2D(const concepts::SpaceOnCells<F>& spc,
             const concepts::Vector<F>& sol,
             const concepts::ElementFormula<F>* a = 0):
    ImplicitResidual<F>(spc, sol), a_(a), sol_(sol){
 
 
 
  }
 
  //starts the error computation
  void compute(){
    computeError_();
  }
 
 
  void addRhs(concepts::LinearForm<F>& lform, Real w = 1.0){
    lforms_.push_back(&lform);
    lws_.push_back(w);
    //DEBUGL(1, "lforms updated : " << lforms_);
  }
 
  void addLhs(concepts::BilinearForm<F>& bform, Real w = 1.0){
    bforms_.push_back(&bform);
    ws_.push_back(w);
    //DEBUGL(1, "bforms updated : " << bforms_);
  }
 
 
  //add information on how to compute the norm of local problem,
  //e.g. energynorm, h1norm
 
  // if no weight is given, it is assumed to be constant  1
 
  //up to now everything in L2 norm
  //
  // example 1)
  // func = Grad, w = 0 , s = 1.0
  // implies norm contribution
  //      || Grad[u_{hp}] ||_{L^2(K)}
  // on each element K.
  //
  // example 2)
  // func = Value,   w = ParsedFormula("(1/sqrt(x^2+y^2))")  s = lambda
  // implies the contribution
  //   lambda * || w * u_{hp} ||_{L^2(K)}
  // on each element K where u_h is the approximate solution
  void setErrorNorm(concepts::ElementFunction<F>& func, const concepts::ElementFormula<F>* w = 0, Real s = 1.0){
    norm_fs_.push_back(&func);
    norm_ws_.push_back(w);
    norm_s_.push_back(s);
 
    if( dynamic_cast<hp2D::Value<F>* >(&func))
      norm_d_.push_back(1);
    else if( dynamic_cast<hp2D::Grad<F>* >(&func) )
      norm_d_.push_back(2);
    else{
      //easy extendable for Curl in 2d add dim = 2, for laplace 1 etc...
      throw conceptsException(concepts::MissingFeature("Other function types currently not setted"));
    }
 
    //DEBUGL(1, "functions updated : "<< norm_fs_)
  }
 
 
 
 
 
  virtual const Real operator()() const {return this->operator()();}
 
  virtual ImplicitResidual2D<F>* clone() const{return new ImplicitResidual2D<F>(*this);}
 
 
protected:
  virtual std::ostream& info(std::ostream& os) const {
    return os << "Implicit Residual Error Estimator 2D()";
  }
 
private:
 
  //isotropic scalar coefficient function in higest derivative operator
  const concepts::ElementFormula<F>* a_;
 
  //the needed information to build local problems, i.e. the lfhs consisting of billinearforms
  concepts::Sequence<concepts::BilinearForm<F>* > bforms_;
  concepts::Sequence<Real> ws_;
 
 
  concepts::Sequence<concepts::LinearForm<F>* > lforms_;
  concepts::Sequence<Real> lws_;
 
 
 
 
//  //norm construction stuff
  concepts::Sequence<concepts::ElementFunction<F>* > norm_fs_;
  concepts::Sequence<const concepts::ElementFormula<F>* > norm_ws_;
  concepts::Sequence<Real> norm_s_;
  concepts::Sequence<uint> norm_d_;
 
 
 
 
 
 
 
  const concepts::Vector<F>& sol_;
 
  virtual void computeError_(){
 
    //mean should be multiplied, thats all from inside
    //the billinearforms must be builded from outside anyway
    if(a_)
      throw conceptsException(concepts::MissingFeature("Coefficient not supported"));
 
 
    //set up needed size fore errors
    LocalEstimator<F>::locError_.resize(this->spc_.nelm());
    concepts::Set<uint> zero; zero.insert(0);
 
 
    //at this time we assume all inner edges have attribute 0.
    hp2D::NeumannTraceSpace ntspc(this->spc_, zero , concepts::EdgeTraceType::MEAN);
    concepts::ElementFormulaVector<1> apprxMean(ntspc,this->sol_, hp2D::Value<Real>());
    //l form due to the apprx Mean
    hp1D::Riesz<Real> lform_mean(apprxMean);
 
 
 
    const hp2D::hpAdaptiveSpaceH1* spc =dynamic_cast<const hp2D::hpAdaptiveSpaceH1*> (&(this->spc_));
    conceptsAssert(spc, concepts::Assertion());
 
    const concepts::BoundaryConditions& bc = *(spc->helper().bc());
 
    DEBUGL(COMPUTE_D, "start iteration");
 
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
    while (!sc->eos()) {
      concepts::ElementWithCell<Real>& element = (*sc)++;
      const hp2D::Quad<Real>* quad = dynamic_cast<const hp2D::Quad<Real>*> (&element);
      conceptsAssert(quad, concepts::Assertion());
 
      concepts::MutableMesh2 msh;
      msh.addCell(&quad->cell(),false);
 
      uint K = quad->support().key();
 
 
      concepts::BoundaryConditions bc_K;
 
      //flag that marks cell intersects a boundary part that is free in sense
      //of dof, i.e no Dirichlet Boundary, but for example a Neumann Boundary
      //bool bnd_free = false; // <=> below set is non empty
      //attributes of edges intersection those boundaries
      concepts::Set<uint> bnd_free_attrb;
 
      const concepts::Edge* edge = 0;
      for(uint i = 0 ; (edge = quad->support().edge(i)) != 0; ++i){
 
        //edge has a non default attribute, e.g. a setted boundary edge
        // that is free in sense of dof, or a setted inner edge
        if(edge->attrib().attrib() != 0){
          //edge is free, so no dirichlet
          if(bc(edge->attrib()).type() != concepts::Boundary::DIRICHLET)
            bnd_free_attrb.insert(edge->attrib().attrib());
          //copy global boundary conditions to the local
          bc_K.add(edge->attrib().attrib(), bc( edge->attrib()).type());
          //FIXME: if inner edges are not zero attributed,
          // on the local problem this inner edge, as boundary edge
          // is marked as FREE, this is ok in the freedom of dofs
          // on that edge, but may lack in consistency
        }
        //edge is a default inner edge
        else{
          //we add neumann condition, as local problems will be defined as
          bc_K.add(edge->attrib().attrib(), concepts::Boundary::NEUMANN);
        }
      }//for edges
 
      //the local polynomial degree, by default + 2
      uint lP = 3;
      const ushort* p = spc->prebuild().innerDof(quad->support());
      if(p){
        conceptsAssert( p[0] == p[1] , concepts::Assertion());
        lP = p[0] + lP - 1 ;
      }
 
      DEBUGL(COMPUTE_D, "on Cell K = " << K << " +++++++++++");
 
      //all neumann conditions
      hp2D::hpAdaptiveSpaceH1 lSpc(msh, 0, lP, &bc_K);
      lSpc.rebuild();
 
      //the current solution vector
      concepts::Vector<Real> sol(lSpc.dim());
      //the rhs vector
      concepts::Vector<Real> rhs(lSpc.dim());
 
 
      DEBUGL(COMPUTE_D, "start build system matrix ");
 
      concepts::SparseMatrix<F> A(lSpc.dim());
      //build local system matrix
      for(uint k = 0 ; k < bforms_.size(); ++k){
        concepts::SparseMatrix<F> A_k(lSpc,*bforms_[k]);
        A_k.addInto(A, ws_[k]);
      }
 
      DEBUGL(COMPUTE_D, "done ...");
 
      //build rhs due to weak residual
      for(uint l = 0 ; l < lforms_.size(); ++l)
        rhs +=  (concepts::Vector<F>(lSpc, *(lforms_[l])) * lws_[l] );
 
      DEBUGL(COMPUTE_D, "rhs weak residual = " << rhs);
 
      //complete rhs with integral about the neumann part.
 
      if(ntspc.nelm()){
      //tracespace on bdn edges, with attribute 0;
      hp2D::TraceSpace tspc_mean(lSpc, zero);
      DEBUGL(COMPUTE_D, " mean tracespace = " << tspc_mean);
      concepts::Vector<F> rhs_mean(tspc_mean.dim());
      std::auto_ptr<hp2D::TraceSpace::Scan> tsc(tspc_mean.scan());
      while ( *(tsc.get()) ) {
 
          hp1D::BaseElement<F>& elm =  (*tsc)++;
 
          DEBUGL(COMPUTE_D, " scan " << elm);
          //hp space and therefore traceSpace are Real
          const concepts::TMatrixBase<Real>& T = elm.T();
          if (T.m() > 0) {
            DEBUGL(COMPUTE_D, " T.m() > 0 ");
            //key of first underlying element
            //DEBUGL(COMPUTE_D, " uelm : " << ntspc.uelm(elm.support()));
 
            bool isFirst = ntspc.isFirstUelm(elm.support(), K);
            DEBUGL(COMPUTE_D, "is First = " << isFirst);
 
            //uint K1 = ntspc.uelm(elm.support())[0].elm->support().key();
            //scalar for the mean, depends on the first underlying element
            Real s = (isFirst) ? 1.0 : -1.0;
            DEBUGL(COMPUTE_D, " s = " << s);
            concepts::ElementMatrix<F> A(1, 1);
            concepts::ElementMatrix<F> B(1, 1);
            DEBUGL(COMPUTE_D, " compute mean application");
            lform_mean(elm, A);
            DEBUGL(COMPUTE_D, " transpose and write into B");
            A.transpose(); T(A, B); B.transpose();
            DEBUGL(COMPUTE_D, " done ...");
            DEBUGL(COMPUTE_D, " B = " << B);
            for(int i = T.n(); i--;){
              DEBUGL(COMPUTE_D, "write into rhs at " << T.index(i) << "  while size = "<< rhs_mean.size());
              rhs_mean[T.index(i)] += s * B(i,0);
            }
          }
          DEBUGL(COMPUTE_D, "current contribution done ...");
      }
 
 
      DEBUGL(COMPUTE_D, "rhs mean residual = " << rhs_mean);
      //update rhs;
      rhs += rhs_mean;
      }//if neumanntrace elements
 
      //Cell boundary intersects Domain boundary, e.g. dirichlet/neumann Boundary
      //or a edge is inner edge, but setted with non-zero attribute
      if(bnd_free_attrb.size()){
        //neumann boundary data or sth else must be added
        conceptsAssert(this->g_.size(), concepts::Assertion());
 
        for(uint i = 0 ; i < this->gAttrbs_.size(); ++i){
            //cell has an edge for which neumann contribution must be taken into account
          if((bnd_free_attrb&&this->gAttrbs_[i]).size()){
 
            //DEBUGL(1, "Found Neumann Cell, compute on it");
 
            hp2D::TraceSpace tspc_g(lSpc, this->gAttrbs_[i]);
            hp1D::Riesz<F> lform_g(*(this->g_[i]));
            //add local neumann data contribution
            rhs += concepts::Vector<F>(tspc_g, lform_g);
          }
        }
      }
 
      DEBUGL(COMPUTE_D, "final rhs  = " << rhs);
 
 
      //solve the problem
      try {
        // Solver
        std::auto_ptr<concepts::Operator<F> > solver(0);
        concepts::ResourceMonitor timer;
        solver.reset(new concepts::Mumps<F>(A));
        conceptsAssert(solver.get(), concepts::Assertion());
        (*solver)(rhs, sol);
      } catch (concepts::ExceptionBase& e) {
        std::cout << e << std::endl;
        throw conceptsException(concepts::MissingFeature("Local problem seems not to be solveable, not supported"));
        return;
      }
 
      //DEBUGL(1, "sol = " << sol);
 
      //set the error of the solution in the requested norm
      Real error_K = 0;
      for(uint n = 0; n < norm_fs_.size(); ++n){
        const uint dim = norm_d_[n];
        if(dim == 1){
          concepts::ElementFormulaVector<1> sol_frm(lSpc, sol, *(norm_fs_[n]));
          if(norm_ws_[n]){
            Real contr = norm_s_[n] * concepts::L2product(lSpc, *(norm_ws_[n]) * sol_frm);
            //DEBUGL(1, " contr(1) = " << contr << "   s = " << norm_s_[n]);
            error_K += contr;
          }
          else
            error_K += norm_s_[n] * concepts::L2product(lSpc, sol_frm);
        }
        if(dim == 2){
          concepts::ElementFormulaVector<2> sol_frm(lSpc, sol, *(norm_fs_[n]));
          if(norm_ws_[n])
            error_K += norm_s_[n] * concepts::L2product(lSpc, *(norm_ws_[n]) * sol_frm);
          else{
            Real contr = norm_s_[n] * concepts::L2product(lSpc, sol_frm);
            //DEBUGL(1, " contr(2) = " << contr);
            error_K += contr;
          }
        }
        //DEBUGL(1, "current error : "<< error_K);
      }
 
      //add current error squared
      this->globError_ += error_K;
 
      //here the global error must be increased and later square rooted
      LocalEstimator<F>::locError_[K] = std::sqrt(error_K);
      DEBUGL(COMPUTE_D, "estimated error on K = " << K << "  : " << LocalEstimator<F>::locError_[K])
 
    }//while
 
    // n^2 = sum n_K^2
    this->globError_ = std::sqrt(this->globError_);
  }
};
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
//namespace hp3D{
//
//class ExplicitResidual3D: concepts::ExplicitResidual{
//public:
//
//  template<typename F>
//    ExplicitResidual3D(const concepts::SpaceOnCells<F>& spc, const concepts::Vector<F>& sol, const concepts::ElementFormula<F>& res);
//
//protected:
//  virtual std::ostream& info(std::ostream& os) const{
//    return os << "Explicit Residual Error Estimator()";
//  }
//
//
//private:
//  virtual void computeError_(const concepts::ElementFormula<F>& res, concepts::HashMap<Real>& jumpResidual);
//  virtual void computeJumpPart_(concepts::HashMap<Real>& jumpResidual) const;
//};
//
//}
 
 
 
 
 
 
} // namespace estimator
 
#endif // implicitResidual_hh