space/integral.hh

Go to the documentation of this file.

/* Integrals of formulas on elements

    Kersten Schmidt, 2005
 */

#ifndef spcIntegral_hh
#define spcIntegral_hh

#include <set>
#include "basics/typedefs.hh"
#include "basics/vectorsMatrices.hh"
#include "basics/vectorsMatricesForward.hh"
#include "geometry/integral.hh"
#include "space/element.hh"
#include "space/formula.hh"
#include "space/postProcess.hh"
#include "space/space.hh"

#define ElemFrmIntegrate_D 0
#define PWintegrate_D      0
#define SpaceIntegrate_D   0
#define ElemFrmL2Product_D 0
#define SpcFrmL2Product_D  0

namespace concepts {

  // forward declaration
  template<typename F, typename G>
  class ElementFormula;


  // ************************************************************* integrate **

  template<typename G>
  Real integrate(const Element<G>& elm)
  {
    const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm);
    if (!cell)
      throw conceptsException(MissingFeature("not a integration cell"));

    Real val(0);
    uint i = 0;
    IntegrationCell::intPoint p;

    while(cell->quadraturePoint(i++, p)) val += p.intermediate;

    return val;
  }

  template<typename F, typename G>
  F integrate(const ElementWithCell<G>& elm, const ElementFormula<F,G>& frm,
              const Real t = 0.0,
              IntegrationCell::intFormType form = IntegrationCell::ZERO) {
    const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm);
    if (!cell)
      throw conceptsException(MissingFeature("not a integration cell"));

    F val(0);
    uint i = 0;
    IntegrationCell::intPoint p;
    DEBUGL(ElemFrmIntegrate_D, "elm = " << elm);
    
    while(cell->quadraturePoint(i++, p, form, true)) {
      DEBUGL(ElemFrmIntegrate_D, "coord = " << p.coord << ", frm(coord) = " <<
             frm(elm, p.coord, t) << ", interme = " << p.intermediate);
      val += frm(elm, p.coord, t) * p.intermediate;
    }
    return val;
  }

  template<class F, typename G>
  F integrate(SpaceOnCells<G>& spc, const ElementFormula<F,G>& frm,
              const Real t = 0.0,
              IntegrationCell::intFormType form = IntegrationCell::ZERO) {
    F val(0);
#if SpaceIntegrate_D
    uint i = 0;
#endif
    // loop over elements
    Scan<ElementWithCell<G> >* sc(spc.scan());
    while (*sc) {
      const ElementWithCell<G>& elm = (*sc)++;
#if SpaceIntegrate_D
      DEBUGL(1, ++i << ".Element: " << elm);
      F tmp = integrate(elm, frm, t, form);
      val += tmp;
      DEBUGL(1, "    Integral: " << tmp);
#else
      val += integrate(elm, frm, t, form);
#endif
    }
    return val;
  }

  // ************************************************************* L2product **

  template<typename F, typename G>
  Real L2product(const ElementWithCell<G>& elm, const ElementFormula<F,G>& u,
                 const ElementFormula<Real>* c = 0, const Real t = 0.0,
                 IntegrationCell::intFormType form = IntegrationCell::ZERO) {
    const IntegrationCell* cell = dynamic_cast<const IntegrationCell*>(&elm);
    if (!cell)
      throw conceptsException(MissingFeature("not a integration cell"));

    DEBUGL(ElemFrmL2Product_D, "Element = " << elm);
    DEBUGL(ElemFrmL2Product_D, "u       = " << u);
    Real val(0);
    uint i = 0;
    IntegrationCell::intPoint p;

    while(cell->quadraturePoint(i++, p, form, true)) {
#if ElemFrmL2Product_D
      if (!c) {
        DEBUGL(1, "coord = " << p.coord << ", u(coord) = " <<
               u(elm, p.coord, t) << ", interme = " << p.intermediate);
      } else {
        DEBUGL(1, "coord = " << p.coord << ", u(coord) = " <<
               u(elm, p.coord, t) << ", interme = " << p.intermediate <<
               ", formula = " << (*c)(elm, p.coord, t) << ", contribution = " <<
               u(elm, p.coord, t)*u(elm, p.coord, t)*
               p.intermediate*(*c)(elm, p.coord, t));
      }
#endif
      F f = u(elm, p.coord, t);
      Real add = std::norm(f) * p.intermediate;
      if (c) add *= (*c)(elm, p.coord, t);
      val += add;
    }
    // one could introduce special handling for piecewise constant
    // formulas which might be zero on some elements or formulas whose
    // support are just some elements

    DEBUGL(ElemFrmL2Product_D, "val = " << val);
    return val;
  }

  template<class F, typename G>
  Real L2product(SpaceOnCells<F>& spc, const G& u,
                 const ElementFormula<Real>* c = 0, const Real t = 0.0,
                 IntegrationCell::intFormType form = IntegrationCell::ZERO) {
    Real val(0);
    DEBUGL(SpcFrmL2Product_D, "spc = " << spc);
    DEBUGL(SpcFrmL2Product_D, "u = " << u);
    // loop over elements
    Scan<ElementWithCell<F> >* sc(spc.scan());
    uint j = 0;
    while (*sc) {
      const ElementWithCell<F>& elm = (*sc)++;
      DEBUGL(SpcFrmL2Product_D, ++j << "th element = " << elm.cell());
      val += L2product(elm, u, c, t, form);
      DEBUGL(SpcFrmL2Product_D, "val = " << val);
    }
    DEBUGL(SpcFrmL2Product_D, "done.");
    return val;
  }

  // ********************************************************** CellIntegral **

  template<class F>
    class CellIntegral : public CellPostprocess<Real> {
      public:
        CellIntegral(F& val, const Formula<F>& formula,
            const std::set<uint>* attributes = 0);
      protected:
        virtual std::ostream& info(std::ostream& os) const;
        F& val_;
        const concepts::Formula<F>& formula_;
        const std::set<uint>* attributes_;
    };

  // ****************************************************** CellFaceIntegral **

  template<class F>
    class CellFaceIntegral : public CellIntegral<F> {
      public:
        CellFaceIntegral(F& val, const concepts::Formula<F>& formula,
            const std::set<uint>* attributes = 0);
        virtual void operator() (const Element<Real>& elm);
        virtual void operator() (const Cell& cell);

      protected:
        virtual std::ostream& info(std::ostream& os) const;
      private:
        template<class G>
          F compute_(const G& elm);
    };

  // ****************************************************** CellEdgeIntegral **

  template<class F>
    class CellEdgeIntegral : public CellIntegral<F> {
      public:
        CellEdgeIntegral(F& val, const concepts::Formula<F>& formula,
            const std::set<uint>* attributes = 0);
        virtual void operator() (const Element<Real>& elm);
        virtual void operator() (const Cell& cell);
      protected:
        virtual std::ostream& info(std::ostream& os) const;
      private:
        template<class G>
          F compute_(const G& elm);
    };

} // namespace concepts

#endif // spcIntegral_hh