app-kersten/ThinSheet/DiffReactAsympModel.hh

Go to the documentation of this file.

#ifndef DiffReactAsympModel_hh
#define DiffReactAsympModel_hh

#include "basics.hh"
#include "toolbox.hh"
#include "geometry.hh"
#include "space.hh"
#include "function.hh"
#include "operator/denseMatrix.hh"
#include "operator/sparseMatrix.hh"
#include "models/models.hh"
#include "hp1D/element.hh"
#include "hp1D/pSpaceL2.hh"
#include "hp2D.hh"

namespace concepts {

  // forward declaration
  template<uint order, class F>
  class DiffReactAsympModelOffset;


  // ******************************************** DiffReactAsympModelSolBase **

  template<class F = Real>
  class DiffReactAsympModelSolBase : public OutputOperator {
  public:
    DiffReactAsympModelSolBase(const F c0);
    ~DiffReactAsympModelSolBase();

    const ElementFormula<F>* operator()(const uint i) const;

    const ElementFormula<Point<F,2> >* grad(const uint i) const;

    const ElementFormula<F>* offset(const uint i) const;

    const ElementFormula<Point<F,2> >* offset_grad(const uint i) const;

    const ElementFormula<F>* source(const uint i) const;


    ElementFormula<F>* trace(const uint i) const;

    ElementFormula<F>* traceDeriv(const uint i) const;

    ElementFormula<F>* meandnu(const uint i) const;

    ElementFormula<F>* U(const uint i) const;

    ElementFormula<Point<F,2> >* grad_U(const uint i) const;

    Space<Real>* space() const {
      conceptsAssert(spc_, Assertion());
      return spc_;
    }

    hp2D::TraceSpace* traceSpace() const {
      conceptsAssert(tspc_, Assertion());
      return tspc_;
    }
    void space(Space<Real>& spc) { spc_ = &spc; }
    void traceSpace(hp2D::TraceSpace& tspc) { tspc_ = &tspc; }

    F c0() const { return c0_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    const F c0_;
    // Offset function from boundary condition
    typename std::auto_ptr<ElementFormula<Real> > u_bc_;
    Space<Real>* spc_;
    hp2D::TraceSpace* tspc_;
    // Sequence of FEM solution vectors uhat0, uhat1, ... in 2D space
    typename concepts::Sequence<Vector<F>*> uvec_;
    // Sequence of FEM solution vectors uhat0, uhat1, ... on midline edges
    typename concepts::Sequence<Vector<F>*> utvec_;
    // Sequence of offset functions utilde0, utilde1, ...
    typename concepts::Sequence<ElementFormula<Real>*> utilde_;
    // Sequence of gradient of the offset function utilde0, utilde1, ...
    typename concepts::Sequence<ElementFormula<Real2d>*> utilde_grad_;
    // The source
    std::auto_ptr<ElementFormula<F> > source_;
    // Sequence of FEM solutions uhat0, uhat1, ...
    typename concepts::Sequence<ElementFormula<F>*> uhat_;
    // Sequence of gradient of FEM solutions uhat0, uhat1, ...
    typename concepts::Sequence<ElementFormula<Point<F,2> >*> uhat_grad_;
    // Sequence of solution as formula vectors u0, u1, ...
    typename concepts::Sequence<ElementFormula<F>*> u_;
    // Sequence of gradient of solution as formula vectors u0, u1, ...
    typename concepts::Sequence<ElementFormula<Point<F,2> >*> u_grad_;
    // Trace of solution
    typename concepts::Sequence<ElementFormula<F>*> ut_;
    // First (tangential) derivative of trace of solution
    typename concepts::Sequence<ElementFormula<F>*> ut_t_;
    // Neumann trace of solution (mean value)
    typename concepts::Sequence<ElementFormula<F>*> meandnu_;
    // Solution inside the sheet in (t,S) coordinates as element formula
    // over the edges 
    typename concepts::Sequence<ElementFormula<F>*> U_;
    // Gradient in (t,S) of the solution inside the sheet in (t,S)
    // coordinates as element formula over the edges
    typename concepts::Sequence<ElementFormula<Point<F,2> >*> grad_U_;
  };

  // ************************************************ DiffReactAsympModelSol **

  template<class F = Real>
  class DiffReactAsympModelSol : public DiffReactAsympModelSolBase<F> {
    friend class DiffReactAsympModelOffset<0,F>;
    friend class DiffReactAsympModelOffset<1,F>;
    friend class DiffReactAsympModelOffset<2,F>;
  public:
    DiffReactAsympModelSol(const F c0);
    DiffReactAsympModelSol(const F c0, const BoundaryConditions bc);
    DiffReactAsympModelSol(const F c0,
                           const ElementFormula<F>& source);

    void set(const Vector<F>& u);
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // ******************************************** DiffReactAsympCollModelSol **

  template<class F = Real>
  class DiffReactAsympCollModelSolBase : public DiffReactAsympModelSolBase<F> {
  public:
    DiffReactAsympCollModelSolBase(const F c0, const Real eps,
                                   const BoundaryConditions* bc = 0);

    inline const Real eps() const { return eps_; }

    inline ElementFormula<F>* meandnuDeriv() const {
      conceptsAssert(meandnuDeriv_.get(), Assertion());
      return meandnuDeriv_.get();
    }

    virtual void set(const Vector<F>& u) = 0;

    Space<Real>* dnu_Space() const {
      conceptsAssert(dnu_spc_, Assertion());
      return dnu_spc_;
    }
    void dnu_Space(Space<Real>& dnu_spc) { dnu_spc_ = &dnu_spc; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    const Real eps_;
    std::auto_ptr<ElementFormula<F> > meandnuDeriv_;
    Space<Real>* dnu_spc_;
  };

  // *************************************** DiffReactAsympCollModelSol_Cont **

  template<uint order, uint alpha = 1, class F = Real>
  class DiffReactAsympCollModelSol_Cont 
    : public DiffReactAsympCollModelSolBase<F> {
  public:
    DiffReactAsympCollModelSol_Cont(const F c0, const Real eps)
      : DiffReactAsympCollModelSolBase<F>(c0, eps) {}
    DiffReactAsympCollModelSol_Cont(const F c0, const Real eps,
                                    const BoundaryConditions bc)
      : DiffReactAsympCollModelSolBase<F>(c0, eps, &bc) {}

    virtual void set(const Vector<F>& u);
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // ************************************ DiffReactAsympCollModelSol_Discont **

  template<uint order, uint alpha = 1, class F = Real>
  class DiffReactAsympCollModelSol_Discont 
    : public DiffReactAsympCollModelSolBase<F> {
  public:
    DiffReactAsympCollModelSol_Discont(const F c0, const Real eps)
      : DiffReactAsympCollModelSolBase<F>(c0, eps) {}
    DiffReactAsympCollModelSol_Discont(const F c0, const Real eps,
                                       const BoundaryConditions bc)
      : DiffReactAsympCollModelSolBase<F>(c0, eps, &bc) {}

    virtual ~DiffReactAsympCollModelSol_Discont();

    virtual void set(const Vector<F>& u);
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::auto_ptr<Vector<F> > meandnu_vec_;
  };


  // ********************************************* DiffReactAsympModelOffset **
  
  template<uint order, class F = Real>
  class DiffReactAsympModelOffset : public ElementFormula<F> {
  public:
    DiffReactAsympModelOffset(DiffReactAsympModelSol<F>& sol);
    virtual ~DiffReactAsympModelOffset() {}
  
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
        const Real2d& p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
        const Real3d& p, const Real t = 0.0) const;
    
    virtual DiffReactAsympModelOffset<order,F>* clone() const { 
      return new DiffReactAsympModelOffset<order,F>
  (sol_, E_o_.get(), E_u_.get(), E_ko_.get(), E_ku_.get());
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympModelSol<F>& sol_;

    mutable std::auto_ptr<const hp2D::H1Extension<Real> >
    E_o_, E_u_, E_ko_, E_ku_;

    DiffReactAsympModelOffset
    (DiffReactAsympModelSol<F>& sol,
     const hp2D::H1Extension<Real>* E_o, const hp2D::H1Extension<Real>* E_u,
     const hp2D::H1Extension<Real>* E_ko, const hp2D::H1Extension<Real>* E_ku) :
      sol_(sol), E_o_(E_o ? E_o->clone() : 0), E_u_(E_u ? E_u->clone() : 0)
      , E_ko_(E_ko ? E_ko->clone() : 0), E_ku_(E_ku ? E_ku->clone() : 0) {}
  };

  // ********************************************* DiffReactAsympModelVector **
  
  template<uint order, class F = Real>
  class DiffReactAsympModelVector : public Vector<F> {
  public:
    DiffReactAsympModelVector(Space<Real>& spc,
                              DiffReactAsympModelSol<F>& sol);
  private:
    const Space<Real>& spc_;
  };

  // *************************************** DiffReactAsympModelSimpleVector **
  
  template<class F = Real>
  class DiffReactAsympModelSimpleVector : public Vector<F> {
  public:
    DiffReactAsympModelSimpleVector(Space<Real>& spc,
                                    DiffReactAsympModelSolBase<F>& sol);
  private:
    const Space<Real>& spc_;
  };

  // ********************************************** DiffReactAsympModelUBase **

  template<typename F = Real>
  class DiffReactAsympModelUBase : public ElementFormula<F> {
  public:
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real2d& p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real3d& p, const Real t = 0.0) const;
    
    virtual Array<F> polynomial
    (const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const = 0;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // ************************************************** DiffReactAsympModelU **

  template<uint order, typename F = Real>
  class DiffReactAsympModelU : public DiffReactAsympModelUBase<F> {
  public:
    DiffReactAsympModelU(DiffReactAsympModelSol<F>& sol)
      : sol_(sol), dtu0_(0), dttu0_(0) {}
    virtual ~DiffReactAsympModelU() {}
      
    virtual DiffReactAsympModelU<order, F>* clone() const { 
      return new DiffReactAsympModelU<order, F>(sol_);
    }

    virtual Array<F> polynomial
    (const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympModelSol<F>& sol_;

    mutable std::auto_ptr<const ElementFormula<F> > dtu0_, dttu0_;
  };

  // ********************************** DiffReactAsympCollModel_U_alpha0and1 **

  template<uint order, uint alpha = 1, typename F = Real>
  class DiffReactAsympCollModelU_alpha0and1 
    : public DiffReactAsympModelUBase<F> {
  public:
    DiffReactAsympCollModelU_alpha0and1
    (DiffReactAsympCollModelSolBase<F>& sol) : sol_(sol) {}
    virtual ~DiffReactAsympCollModelU_alpha0and1() {}
      
    virtual DiffReactAsympCollModelU_alpha0and1<order,alpha,F>* clone() const {
      return new DiffReactAsympCollModelU_alpha0and1<order,alpha,F>(sol_);
    }

    virtual Array<F> polynomial
    (const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympCollModelSolBase<F>& sol_;
  };

  // ************************************** DiffReactAsympCollModel_U_alpha2 **

  template<uint order, typename F = Real>
  class DiffReactAsympCollModelU_alpha2
    : public DiffReactAsympModelUBase<F> {
  public:
    DiffReactAsympCollModelU_alpha2
    (DiffReactAsympCollModelSolBase<F>& sol) : sol_(sol) {}
    virtual ~DiffReactAsympCollModelU_alpha2() {}
      
    virtual DiffReactAsympCollModelU_alpha2<order,F>* clone() const {
      return new DiffReactAsympCollModelU_alpha2<order,F>(sol_);
    }

    F operator() (const ElementWithCell<Real>& elm,
                  const Real2d& p, const Real t = 0.0) const;

    virtual Array<F> polynomial
    (const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympCollModelSolBase<F>& sol_;
  };

  // ********************************************** DiffReactAsympCollModelU **

  template<uint order, uint alpha = 1, typename F = Real>
  class DiffReactAsympCollModelU {};

  template<uint order, typename F>
  class DiffReactAsympCollModelU<order,0,F> : 
    public DiffReactAsympCollModelU_alpha0and1<order,0,F> {  
  public:
    DiffReactAsympCollModelU
    (DiffReactAsympCollModelSolBase<F>& sol)
      : DiffReactAsympCollModelU_alpha0and1<order,0,F>(sol) {}
  };

  template<uint order, typename F>
  class DiffReactAsympCollModelU<order,1,F> : 
    public DiffReactAsympCollModelU_alpha0and1<order,1,F> {  
  public:
    DiffReactAsympCollModelU
    (DiffReactAsympCollModelSolBase<F>& sol) 
      : DiffReactAsympCollModelU_alpha0and1<order,1,F>(sol) {}
  };

  template<uint order, typename F>
  class DiffReactAsympCollModelU<order,2,F> : 
    public DiffReactAsympCollModelU_alpha2<order,F> {  
  public:
    DiffReactAsympCollModelU
    (DiffReactAsympCollModelSolBase<F>& sol) 
      : DiffReactAsympCollModelU_alpha2<order,F>(sol) {}
  };

  // ********************************************** DiffReactAsympModelUepsN **

  template<typename F = Real>
  class DiffReactAsympModelUepsN : public ElementFormula<F> {
  public:
    DiffReactAsympModelUepsN(DiffReactAsympModelSol<F>& sol, const Real eps,
                             uint order) 
      : sol_(sol), eps_(eps), order_(order) {}
    virtual ~DiffReactAsympModelUepsN() {}
  
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real2d& p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real3d& p, const Real t = 0.0) const;
    
    virtual DiffReactAsympModelUepsN<F>* clone() const { 
      return new DiffReactAsympModelUepsN<F>(sol_, eps_, order_);
    }

    Real eps() const { return eps_; }
    uint order() const { return order_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympModelSol<F>& sol_;
    const Real eps_;
    const uint order_;
  };

  // ********************************************** DiffReactAsympModelUBase **

  template<typename F = Real>
  class DiffReactAsympModelgrad_UBase : public ElementFormula<Point<F,2> > {
  public:
    virtual Point<F,2> operator() (const ElementWithCell<Real>& elm,
                                   const Real p, const Real t = 0.0) const;
    virtual Point<F,2> operator() (const ElementWithCell<Real>& elm,
                                   const Real2d& p, const Real t = 0.0) const;
    virtual Point<F,2> operator() (const ElementWithCell<Real>& elm,
                                   const Real3d& p, const Real t = 0.0) const;
    
    virtual Array<F> polynomial
    (bool normal, const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const = 0;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };

  // ***************************** DiffReactAsympCollModel_grad_U_alpha0and1 **

  template<uint order, uint alpha = 1, typename F = Real>
  class DiffReactAsympCollModel_grad_U_alpha0and1 
    : public DiffReactAsympModelgrad_UBase<F> {
  public:
    DiffReactAsympCollModel_grad_U_alpha0and1
    (DiffReactAsympCollModelSolBase<F>& sol) : sol_(sol) {}
    virtual ~DiffReactAsympCollModel_grad_U_alpha0and1() {}
      
    virtual DiffReactAsympCollModel_grad_U_alpha0and1<order,alpha,F>* clone() const {
      return new DiffReactAsympCollModel_grad_U_alpha0and1<order,alpha,F>(sol_);
    }

    virtual Array<F> polynomial
    (bool normal, const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympCollModelSolBase<F>& sol_;
  };

  // ********************************* DiffReactAsympCollModel_grad_U_alpha2 **

  template<uint order, typename F = Real>
  class DiffReactAsympCollModel_grad_U_alpha2 
    : public DiffReactAsympModelgrad_UBase<F> {
  public:
    DiffReactAsympCollModel_grad_U_alpha2
    (DiffReactAsympCollModelSolBase<F>& sol) : sol_(sol) {}
    virtual ~DiffReactAsympCollModel_grad_U_alpha2() {}
      
    virtual DiffReactAsympCollModel_grad_U_alpha2<order,F>* clone() const {
      return new DiffReactAsympCollModel_grad_U_alpha2<order,F>(sol_);
    }

    Point<F,2> operator() (const ElementWithCell<Real>& elm,
                           const Real2d& p, const Real t = 0.0) const;

    virtual Array<F> polynomial
    (bool normal, const ElementWithCell<Real>& elm,
     const Real p, const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympCollModelSolBase<F>& sol_;
  };

  // ***************************************** DiffReactAsympCollModelgrad_U **

  template<uint order, uint alpha = 1, typename F = Real>
  class DiffReactAsympCollModelgrad_U {};

  template<uint order, typename F>
  class DiffReactAsympCollModelgrad_U<order,0,F> : 
    public DiffReactAsympCollModel_grad_U_alpha0and1<order,0,F> {  
  public:
    DiffReactAsympCollModelgrad_U
    (DiffReactAsympCollModelSolBase<F>& sol)
      : DiffReactAsympCollModel_grad_U_alpha0and1<order,0,F>(sol) {}
  };

  template<uint order, typename F>
  class DiffReactAsympCollModelgrad_U<order,1,F> : 
    public DiffReactAsympCollModel_grad_U_alpha0and1<order,1,F> {  
  public:
    DiffReactAsympCollModelgrad_U
    (DiffReactAsympCollModelSolBase<F>& sol) 
      : DiffReactAsympCollModel_grad_U_alpha0and1<order,1,F>(sol) {}
  };

  template<uint order, typename F>
  class DiffReactAsympCollModelgrad_U<order,2,F> : 
    public DiffReactAsympCollModel_grad_U_alpha2<order,F> {  
  public:
    DiffReactAsympCollModelgrad_U
    (DiffReactAsympCollModelSolBase<F>& sol) 
      : DiffReactAsympCollModel_grad_U_alpha2<order,F>(sol) {}
  };

  // ****************************************** diffReactAsympModelPowerLoss **

  template<class F>
  Real diffReactAsympModelPowerLoss(Real c0, Space<Real>* spc,
                                    DiffReactAsympModelUepsN<F>& UepsN)
  {
    Real val = 0.0;
    hp2D::TraceSpace* tspc = dynamic_cast<hp2D::TraceSpace*>(spc);
    conceptsAssert(tspc, Assertion());
    
    // loop over elements
    std::auto_ptr<hp2D::TraceSpace::Scan> tsc(tspc->scan());
    while (*tsc)
      val += diffReactAsympModelPowerLoss(c0, (*tsc)++, UepsN);
    return val;
  }

  template<class F>
  Real diffReactAsympModelPowerLoss(Real c0, hp1D::Element<Real>& elm,
                                    DiffReactAsympModelUepsN<F>& UepsN) 
  {
    const ElementWithCell<Real>* e = 
      dynamic_cast<const ElementWithCell<Real>*>(&elm);
    conceptsAssert(e, Assertion());

    const Edge2d* edge = dynamic_cast<const Edge2d*>(&e->cell());
    if (!edge)
      throw conceptsException(MissingFeature("not a valid cell"));

    // rule for integration over S
    QuadratureRuleGaussJacobi integrationS(UepsN.order() + 2);
    
    // quadrature points in t and S direction
    const uint        n[2] = { elm.integration()->n(), integrationS.n() };
    const Real* const w[2] = { elm.integration()->weights(),
                                integrationS.weights() };
    const Real* const q[2] = { elm.integration()->abscissas(),
                                integrationS.abscissas() };

    F val = 0.0;
    // loop over integration points along the midline
    for (uint i = 0; i < n[0]; ++i) {
      // integration point
      Real2d P; P[0] = (q[0][i]+1.0)/2.0;
      // curvature times thickness
      //        | sign ??
      Real kd = -edge->map()->curvature(P[0]) * UepsN.eps();
      //                                |
      //                                wrong position ???
      //                                would be only visible at ellipse

      // loop over integration points in S direction
      F tmp = 0.0;
      for(uint j = 0; j < n[1]; ++j) {
        P[1] = q[1][j] * 1 / 2.;
  tmp += std::norm(UepsN(elm, P)) * w[1][j] * (1 + P[1]*kd);
      } // for j
      val += tmp * w[0][i] * edge->jacobianDeterminant(P[0]);
    } // for i
    // factor 4 due to integration rule over [-1,1] instead of
    // [-0.5,0.5] in S direction and [-1,1] instead of [0,1] in
    // tangential direction
    val *= c0 / 4.0; 
    
    return val;
  }

  // ********************************************** DiffReactAsympModeluepsN **

  template<typename F = Real>
  class DiffReactAsympModeluepsN : public ElementFormula<F> {
  public:
    DiffReactAsympModeluepsN(DiffReactAsympModelSolBase<F>& sol, 
                             CellToCellMapping<2,2>& map, 
                             Space<Real>& spc, Space<Real>& tspc,
                             const Real eps,
                             uint order);
    virtual ~DiffReactAsympModeluepsN() {}
  
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real2d& p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real3d& p, const Real t = 0.0) const;
    
    virtual DiffReactAsympModeluepsN<F>* clone() const { 
      return new DiffReactAsympModeluepsN<F>(sol_, map_, spc_, tspc_, eps_,
                                             order_);
    }

    Real eps() const { return eps_; }
    uint order() const { return order_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympModelSolBase<F>& sol_;
    CellToCellMapping<2,2>& map_;
    Space<Real> &spc_, &tspc_;
    const Real eps_;
    const uint order_;
    
    std::map<const Connector*, const ElementWithCell<Real>*> mapToElm_;
  };


  // ***************************************** DiffReactAsympModeluepsN_grad **

  template<typename F = Real>
  class DiffReactAsympModeluepsN_grad : public ElementFormula<Point<F,2> > {
  public:
    DiffReactAsympModeluepsN_grad(DiffReactAsympModelSolBase<F>& sol, 
                                  CellToCellMapping<2,2>& map, 
                                  Space<Real>& spc, Space<Real>& tspc,
                                  const Real eps,
                                  uint order);
    virtual ~DiffReactAsympModeluepsN_grad() {}
  
    virtual Point<F,2> operator() (const ElementWithCell<Real>& elm,
                                   const Real p, const Real t = 0.0) const;
    virtual Point<F,2> operator() (const ElementWithCell<Real>& elm,
                                   const Real2d& p, const Real t = 0.0) const;
    virtual Point<F,2> operator() (const ElementWithCell<Real>& elm,
                                   const Real3d& p, const Real t = 0.0) const;
    
    virtual DiffReactAsympModeluepsN_grad<F>* clone() const { 
      return new DiffReactAsympModeluepsN_grad<F>(sol_, map_, spc_, tspc_, eps_,
                                                  order_);
    }

    Real eps() const { return eps_; }
    uint order() const { return order_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    DiffReactAsympModelSolBase<F>& sol_;
    CellToCellMapping<2,2>& map_;
    Space<Real> &spc_, &tspc_;
    const Real eps_;
    const uint order_;
    
    std::map<const Connector*, const ElementWithCell<Real>*> mapToElm_;
  };


  // ****************************************** ElementFormulaTimesCurvature **

  template<typename F = Real>
  class ElementFormulaTimesCurvature : public ElementFormula<F> {
  public:
    ElementFormulaTimesCurvature(ElementFormula<F>& frm) : frm_(frm.clone()) {}
    virtual ~ElementFormulaTimesCurvature() {}
  
    virtual F operator() (const ElementWithCell<Real>& elm,
                          const Real p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
        const Real2d& p, const Real t = 0.0) const;
    virtual F operator() (const ElementWithCell<Real>& elm,
        const Real3d& p, const Real t = 0.0) const;
    
    virtual ElementFormulaTimesCurvature<F>* clone() const { 
      return new ElementFormulaTimesCurvature<F>(*frm_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::auto_ptr<ElementFormula<F> > frm_;
  };

  // ********************************************* DiffReactAsympCollFormula **
  
  template<class F>
  class DiffReactAsympCollFormula : public ElementFormula<F> {
  public:
    DiffReactAsympCollFormula(const F c0, const Real eps)
      : c0_(c0), eps_(eps) {}
    F&    c0 () { return c0_; }
    Real& eps() { return eps_;}
    virtual F operator()(const ElementWithCell<Real>& elm, const Real p,
                         const Real t = 0.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;

    virtual F compute_(const Edge2d& cell, const Real p) const = 0;
    
    virtual std::string str_() const = 0;

    F c0_;
    Real eps_;
  };

  // ********************************************* DiffReactAsympCollGamma_0 **

  template<class F, uint N, uint d = 0>
  class DiffReactAsympCollGamma_0 : public DiffReactAsympCollFormula<F> {
  public:
    DiffReactAsympCollGamma_0(const F c0, const Real eps)
      : DiffReactAsympCollFormula<F>(c0, eps) {}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    virtual F compute_(const Edge2d& cell, const Real p) const;

    virtual std::string str_() const = 0;
  };

  // ********************************************* DiffReactAsympCollGamma_1 **

  template<class F, uint N, uint d = 0>
  class DiffReactAsympCollGamma_1 : public DiffReactAsympCollFormula<F> {
  public:
    DiffReactAsympCollGamma_1(const F c0, const Real eps)
      : DiffReactAsympCollFormula<F>(c0, eps) {}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    virtual F compute_(const Edge2d& cell, const Real p) const;

    virtual std::string str_() const = 0;
  };

  // ********************************************* DiffReactAsympCollDelta_0 **

  template<class F, uint N, uint d = 0>
  class DiffReactAsympCollDelta_0 : public DiffReactAsympCollFormula<F> {
  public:
    DiffReactAsympCollDelta_0(const F c0, const Real eps)
      : DiffReactAsympCollFormula<F>(c0, eps) {}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    virtual F compute_(const Edge2d& cell, const Real p) const;

    virtual std::string str_() const = 0;
  };

  // ********************************************* DiffReactAsympCollDelta_1 **

  template<class F, uint N, uint d = 0>
  class DiffReactAsympCollDelta_1 : public DiffReactAsympCollFormula<F> {
  public:
    DiffReactAsympCollDelta_1(const F c0, const Real eps)
      : DiffReactAsympCollFormula<F>(c0, eps) {}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    virtual F compute_(const Edge2d& cell, const Real p) const;

    virtual std::string str_() const = 0;
  };

  // *********************************************** DiffReactAsympModelBase **

  template<class F = Real>
  class DiffReactAsympModelBase : public Model<F> {
  public:
    DiffReactAsympModelBase(std::string problemName) 
      : Model<F>(problemName), c_(1.0), eps_(1.0) {}

    virtual DiffReactAsympModelSolBase<F>* modelSol() = 0;

    inline const Real eps() const { return eps_; }
    inline void setEps(Real eps) {
      if (eps != eps_) {
        this->sol_.reset(0);
        eps_ = eps;
      }
    }
    inline const F& c() const { return c_; }
    inline void setC(const F c) { 
      if (c != c_) {
        this->sol_.reset(0);
        c_ = c;
      }
    }
  protected:
    F c_;
    Real eps_;
  };

  // *********************************************** DiffReactAsympCollModel **

  template<class F = Real, uint order = 0, uint alpha = 1,
           class G = typename Realtype<F>::type>
  class DiffReactAsympCollModel {
  private:
    
  };

  // ***************************** DiffReactAsympCollModel_Alpha0_Order0and1 **

  template<class F = Real, uint order = 2, 
           class G = typename Realtype<F>::type>
  class DiffReactAsympCollModel_Alpha0_Order0and1 
    : public DiffReactAsympModelBase<F> {
  public:
    DiffReactAsympCollModel_Alpha0_Order0and1
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm,
     BoundaryConditions* bc = 0);

    virtual ~DiffReactAsympCollModel_Alpha0_Order0and1();

    virtual Space<Real>& space() const {
      conceptsAssert(spcE_.get(), Assertion());
      return *spcE_;
    }

    virtual DiffReactAsympCollModelSol_Cont<order,0,F>* modelSol() {
      return modelSol_.get();
    }
  private:
    Mesh2& msh_;
    // Space for E-field
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1> spcE_;
    // Trace space on r = Ri
    std::auto_ptr<hp2D::TraceSpace> tspc_;

    std::auto_ptr<DiffReactAsympCollModelSol_Cont<order,0,F> > modelSol_;

    std::auto_ptr<SparseMatrix<Real> > A_;
    std::auto_ptr<SparseMatrix<Real> > M_;
    std::auto_ptr<DenseMatrix<Real> > L_;

    std::auto_ptr<Vector<G> > rhs_;

    virtual Mesh& mesh_() {
      return msh_;
    }

    virtual const std::string mshAbbr_() {
      return "ThinSheet";
    }

    virtual void solve_();
  };


  // ***************************** DiffReactAsympCollModel_Alpha1_Order0and1 **

  template<class F = Real, uint order = 2, 
           class G = typename Realtype<F>::type>
  class DiffReactAsympCollModel_Alpha1_Order0and1 
    : public DiffReactAsympModelBase<F> {
  public:
    DiffReactAsympCollModel_Alpha1_Order0and1
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm,
     BoundaryConditions* bc = 0);

    virtual ~DiffReactAsympCollModel_Alpha1_Order0and1();

    virtual Space<Real>& space() const {
      conceptsAssert(spcE_.get(), Assertion());
      return *spcE_;
    }

    virtual DiffReactAsympCollModelSol_Cont<order,1,F>* modelSol() {
      return modelSol_.get();
    }
  private:
    Mesh2& msh_;
    // Space for E-field
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1> spcE_;
    // Trace space on r = Ri
    std::auto_ptr<hp2D::TraceSpace> tspc_;

    std::auto_ptr<DiffReactAsympCollModelSol_Cont<order,1,F> > modelSol_;

    std::auto_ptr<SparseMatrix<Real> > A_;
    std::auto_ptr<SparseMatrix<Real> > M_;

    std::auto_ptr<Vector<G> > rhs_;

    virtual Mesh& mesh_() {
      return msh_;
    }

    virtual const std::string mshAbbr_() {
      return "ThinSheet";
    }

    virtual void solve_();
  };

  // ***************************** DiffReactAsympCollModel_Alpha1_Order2plus **

  template<class F = Real, uint order = 2, 
           class G = typename Realtype<F>::type>
  class DiffReactAsympCollModel_Alpha1_Order2plus 
    : public DiffReactAsympModelBase<F> {
  public:
    DiffReactAsympCollModel_Alpha1_Order2plus
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm,
     BoundaryConditions* bc = 0);

    virtual ~DiffReactAsympCollModel_Alpha1_Order2plus();

    virtual Space<Real>& space() const {
      conceptsAssert(spcE_.get(), Assertion());
      return *spcE_;
    }

    virtual DiffReactAsympCollModelSol_Discont<order,1,F>* modelSol() {
      return modelSol_.get();
    }

    const Vector<F>* solL() {
      this->solve_();
      return solL_.get();
    }
  private:
    Mesh2& msh_;
    MutableMesh1 mshG_;
    // Space for E-field
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1_DG> spcE_;
    // Trace space for mean {e} on r = Ri
    std::auto_ptr<hp2D::TraceSpace> spcE_mean_;
    std::auto_ptr<hp2D::TraceSpace> spcE_jump_;
    // Space for Lagrange multiplier space lambda = {dn e}
    std::auto_ptr<hp1D::pSpaceL2> spcL_;

    std::auto_ptr<DiffReactAsympCollModelSol_Discont<order,1,F> > modelSol_;

    std::auto_ptr<SparseMatrix<Real> > A_;
    std::auto_ptr<SparseMatrix<Real> > M_;
    std::auto_ptr<SparseMatrix<Real> > E_;
    std::auto_ptr<SparseMatrix<Real> > T01_;
    std::auto_ptr<SparseMatrix<Real> > L_;
    std::auto_ptr<SparseMatrix<Real> > N_;
    std::auto_ptr<SparseMatrix<Real> > Mk2_;

    std::auto_ptr<Vector<G> > rhs_;

    std::auto_ptr<Vector<F> > solL_;

    virtual Mesh& mesh_() {
      return msh_;
    }

    virtual const std::string mshAbbr_() {
      return "ThinSheet";
    }

    virtual void solve_();
  };


  // ********************************* DiffReactAsympCollModel_Alpha2_Order0 **

  template<class F = Real, class G = typename Realtype<F>::type>
  class DiffReactAsympCollModel_Alpha2_Order0 
    : public DiffReactAsympModelBase<F> {
  public:
    DiffReactAsympCollModel_Alpha2_Order0
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm,
     BoundaryConditions* bc = 0, bool robuster = false);

    virtual ~DiffReactAsympCollModel_Alpha2_Order0();

    virtual Space<Real>& space() const {
      conceptsAssert(spcE_.get(), Assertion());
      return *spcE_;
    }

    virtual DiffReactAsympCollModelSol_Cont<0,2,F>* modelSol() {
      return modelSol_.get();
    }
  private:
    Mesh2& msh_;
    // Space for E-field
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1> spcE_;
    // Trace space on r = Ri (not used)
    std::auto_ptr<hp2D::TraceSpace> tspc_;

    bool robuster_;

    std::auto_ptr<DiffReactAsympCollModelSol_Cont<0,2,F> > modelSol_;

    std::auto_ptr<SparseMatrix<Real> > A_;

    std::auto_ptr<SparseMatrix<Real> > M_;

    std::auto_ptr<Vector<G> > rhs_;

    virtual Mesh& mesh_() {
      return msh_;
    }

    virtual const std::string mshAbbr_() {
      return "ThinSheet";
    }

    virtual void solve_();
  };


  // ***************************** DiffReactAsympCollModel_Alpha2_Order1plus **

  template<class F = Real, uint order = 1, class G = typename Realtype<F>::type>
  class DiffReactAsympCollModel_Alpha2_Order1plus 
    : public DiffReactAsympModelBase<F> {
  public:
    DiffReactAsympCollModel_Alpha2_Order1plus
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm,
     BoundaryConditions* bc = 0, bool robuster = false);

    virtual ~DiffReactAsympCollModel_Alpha2_Order1plus();

    virtual Space<Real>& space() const {
      conceptsAssert(spcE_.get(), Assertion());
      return *spcE_;
    }
    virtual DiffReactAsympCollModelSol_Discont<order,2,F>* modelSol() {
      return modelSol_.get();
    }
  private:
    Mesh2& msh_;
    MutableMesh1 mshG_;
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1_DG> spcE_;
    std::auto_ptr<hp2D::TraceSpace> spcE_mean_;
    std::auto_ptr<hp2D::TraceSpace> spcE_jump_;
    std::auto_ptr<hp1D::pSpaceL2> spcL_;

    bool robuster_;

    std::auto_ptr<DiffReactAsympCollModelSol_Discont<order,2,F> > modelSol_;

    std::auto_ptr<SparseMatrix<Real> > A_;
    std::auto_ptr<SparseMatrix<Real> > E_;
    std::auto_ptr<SparseMatrix<Real> > T01_;
    std::auto_ptr<SparseMatrix<Real> > T10_;
    std::auto_ptr<SparseMatrix<Real> > T21_;
    std::auto_ptr<SparseMatrix<Real> > L_;

    std::auto_ptr<Vector<G> > rhs_;

    std::auto_ptr<Vector<F> > solL_;

    virtual Mesh& mesh_() {
      return msh_;
    }

    virtual const std::string mshAbbr_() {
      return "ThinSheet";
    }

    virtual void solve_();
  };


  // ************************************** DiffReactAsympCollModel<F,0-1,0> **

  template<class F, class G>
  class DiffReactAsympCollModel<F, 0, 0, G> 
    : public DiffReactAsympCollModel_Alpha0_Order0and1<F,0, G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0)
      : DiffReactAsympCollModel_Alpha0_Order0and1<F,0,G>(msh, l, p, j0frm, bc) {}
  };

  template<class F, class G>
  class DiffReactAsympCollModel<F, 1, 0, G> 
    : public DiffReactAsympCollModel_Alpha0_Order0and1<F,1,G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0)
      : DiffReactAsympCollModel_Alpha0_Order0and1<F,1,G>(msh, l, p, j0frm, bc) {}
  };


  // ************************************** DiffReactAsympCollModel<F,0-3,1> **

  template<class F, class G>
  class DiffReactAsympCollModel<F, 0, 1, G> 
    : public DiffReactAsympCollModel_Alpha1_Order0and1<F,0,G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0)
      : DiffReactAsympCollModel_Alpha1_Order0and1<F,0,G>(msh, l, p, j0frm, bc) {}
  };

  template<class F, class G>
  class DiffReactAsympCollModel<F, 1, 1, G> 
    : public DiffReactAsympCollModel_Alpha1_Order0and1<F,1,G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0)
      : DiffReactAsympCollModel_Alpha1_Order0and1<F,1,G>(msh, l, p, j0frm, bc) {}
  };

  template<class F, class G>
  class DiffReactAsympCollModel<F, 2, 1, G> 
    : public DiffReactAsympCollModel_Alpha1_Order2plus<F,2,G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0)
      : DiffReactAsympCollModel_Alpha1_Order2plus<F,2,G>(msh, l, p, j0frm, bc) {}
  };

  template<class F, class G>
  class DiffReactAsympCollModel<F, 3, 1, G> 
    : public DiffReactAsympCollModel_Alpha1_Order2plus<F,3,G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0)
      : DiffReactAsympCollModel_Alpha1_Order2plus<F,3,G>(msh, l, p, j0frm, bc) {}
  };

  // ************************************** DiffReactAsympCollModel<F,0-1,2> **

  template<class F, class G>
  class DiffReactAsympCollModel<F, 0, 2, G>
    : public DiffReactAsympCollModel_Alpha2_Order0<F, G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0, bool robuster = false)
      : DiffReactAsympCollModel_Alpha2_Order0<F,G>(msh, l, p, j0frm, bc, robuster) {}
  };

  template<class F, class G>
  class DiffReactAsympCollModel<F, 1, 2, G>
    : public DiffReactAsympCollModel_Alpha2_Order1plus<F,1,G> {
  public:
    DiffReactAsympCollModel
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm, 
     BoundaryConditions* bc = 0, bool robuster = false)
      : DiffReactAsympCollModel_Alpha2_Order1plus<F,1,G>(msh, l, p, j0frm, bc, robuster)
    {}
  };


  // *************************************** DiffReactAsympModelKraehenbuehl **

  template<class F = Real, class G = typename Realtype<F>::type>
  class DiffReactAsympModelKraehenbuehl : public DiffReactAsympModelBase<F> {
  public:
    DiffReactAsympModelKraehenbuehl
    (Mesh2& msh, const uint l, const uint p,
     const ElementFormula<G>& j0frm,
     BoundaryConditions* bc = 0);

    virtual ~DiffReactAsympModelKraehenbuehl();

    virtual Space<Real>& space() const {
      conceptsAssert(spcE_.get(), Assertion());
      return *spcE_;
    }
    virtual DiffReactAsympCollModelSol_Discont<0,0,F>* modelSol() {
      return modelSol_.get();
    }
  private:
    Mesh2& msh_;
    MutableMesh1 mshG_;
    std::auto_ptr<hp2D::hpAdaptiveSpaceH1_DG> spcE_;
    std::auto_ptr<hp2D::TraceSpace> spcE_mean_;
    std::auto_ptr<hp2D::TraceSpace> spcE_jump_;
    std::auto_ptr<hp1D::pSpaceL2> spcL_;

    // second template parameter is dummy as no internal solution given
    std::auto_ptr<DiffReactAsympCollModelSol_Discont<0,0,F> > modelSol_;

    std::auto_ptr<SparseMatrix<Real> > A_;
    std::auto_ptr<SparseMatrix<Real> > E_;
    std::auto_ptr<SparseMatrix<Real> > T01_;
    std::auto_ptr<SparseMatrix<Real> > L_;

    std::auto_ptr<Vector<G> > rhs_;

    std::auto_ptr<Vector<F> > solL_;

    virtual Mesh& mesh_() {
      return msh_;
    }

    virtual const std::string mshAbbr_() {
      return "ThinSheet";
    }

    virtual void solve_();
  };

    
} // namespace concepts

#endif // DiffReactAsympModel_hh