elementMaps.hh

Go to the documentation of this file.

 
#ifndef elementMaps_hh
#define elementMaps_hh
 
#include "basics/outputOperator.hh"
#include "basics/typedefs.hh"
#include "basics/vectorsMatrices.hh"
#include "basics/cloneable.hh"
#include "toolbox/sequence.hh"
#include "toolbox/arrayOp.hh"
 
namespace concepts {
 
  // ***************************************************************** Map1d **
 
  class Map1d : public OutputOperator {
  public:
    std::ostream& info(std::ostream& os) const {return os << concepts::typeOf(*this)<<"()";};
  };
 
  // ***************************************************************** Map2d **
 
  class Map2d : public OutputOperator {
  public:
    std::ostream& info(std::ostream& os) const {return os << concepts::typeOf(*this)<<"()";};
    virtual bool straight() const {return false;};
  };
 
  // ************************************************************* MapEdge1d **
 
  class MapEdge1d : public Map1d {
  public:
    inline MapEdge1d(const Real vtx0, const Real vtx1)  :
      vtx0_(vtx0), vtx1_(vtx1) {}
 
    inline Real operator()(const Real x) const {
      return vtx0_*(1.0-x) + vtx1_*x; }
 
    inline Real jacobian(const Real) const { return vtx1_-vtx0_; }
 
    inline MapEdge1d* clone() const {return new MapEdge1d(vtx0_, vtx1_);};
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Real vtx0_, vtx1_;
  };
 
  // ********************************************************* MappingEdge2d **
 
  class MappingEdge2d : public Cloneable, public Map1d {
  public:
    MappingEdge2d(const Real2d vtx0, const Real2d vtx1);
    MappingEdge2d(const MappingEdge2d &edgemap);
    virtual MappingEdge2d* clone() const = 0;
    virtual MappingEdge2d* inverse() const;
    virtual Real2d operator()(const Real t) const = 0;
    virtual Real2d derivative(const Real t, const uint n = 1) const = 0;
 
    virtual Real2d normal(const Real t) const;
    inline Real2d n0(const Real t) const {
      return UnitNd<2>(normal(t));
    }
    virtual Real curvature(const Real t, uint n = 0) const;
 
    virtual MappingEdge2d* part(const Real t0, const Real t1) const;
 
    inline const Real2d& vtx(uint i) const {return vtx_[i];}
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    Real2d vtx_[2];
  };
 
  // ***************************************************** PartMappingEdge2d **
 
  class PartMappingEdge2d : public MappingEdge2d {
  public:
    PartMappingEdge2d(const MappingEdge2d& edgemap, 
                      const Real t0, const Real t1);
    PartMappingEdge2d(const PartMappingEdge2d &edgemap);
    virtual PartMappingEdge2d* clone() const;
    virtual PartMappingEdge2d* inverse() const;
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::unique_ptr<const MappingEdge2d> edgemap_;
    const Real t0_, t1_;
 
    inline const Real t_(const Real t) const {
      return t0_*(1. - t) + t1_*t;
    }
  };
 
  // ************************************************** InverseMappingEdge2d **
 
  class InverseMappingEdge2d : public MappingEdge2d {
  public:
    InverseMappingEdge2d(const MappingEdge2d& edgemap);
 
    virtual InverseMappingEdge2d* clone() const;
    virtual MappingEdge2d* inverse() const {
      return edgemap_->clone();
    }
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    virtual Real curvature(const Real t, uint n = 0) const;
    PartMappingEdge2d* part(const Real t0, const Real t1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    const MappingEdge2d* edgemap_;
  };
 
  // ************************************************* MappingStraightEdge2d **
 
  class MappingStraightEdge2d : public MappingEdge2d {
  public:
    MappingStraightEdge2d(const Real2d vtx0, const Real2d vtx1);
    MappingStraightEdge2d(const MappingStraightEdge2d &edgemap);
    virtual MappingStraightEdge2d* clone() const;
    virtual MappingStraightEdge2d* inverse() const {
      return new MappingStraightEdge2d(vtx(1), vtx(0));
    }
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    virtual Real curvature(const Real t, uint n = 0) const;
    MappingStraightEdge2d* part(const Real t0, const Real t1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    Real2d delta_;
  };
 
  // ************************************************** CircleMappingEdge2d **
 
  class CircleMappingEdge2d : public MappingEdge2d {
  public:
    CircleMappingEdge2d(const Real r, const Real2d vtx0, const Real2d vtx1);
    CircleMappingEdge2d(const Real2d center, const Real r,
                        const Real angle0, const Real angle1);
    CircleMappingEdge2d(const Real2d vtx0, const Real2d vtxm, const Real2d vtx1);
    CircleMappingEdge2d(const CircleMappingEdge2d &edgemap);
    virtual CircleMappingEdge2d* clone() const;
    virtual MappingEdge2d* inverse() const;
 
    bool straight() const;
 
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    virtual Real curvature(const Real t, const uint n = 0) const;
    MappingEdge2d* part(const Real t0, const Real t1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Real r_;
    Real2d m_;
    Real angle0_, angle1_;
 
    inline Real phi_(const Real t) const {
      return angle0_ * (1.-t) + angle1_ * t;
    }
  };
 
  // ************************************************** EllipseMappingEdge2d **
 
  class EllipseMappingEdge2d : public MappingEdge2d {
  public:
    EllipseMappingEdge2d(const Real2d center, const Real a, const Real b,
                         const Real2d vtx0, const Real2d vtx1);
    EllipseMappingEdge2d(const EllipseMappingEdge2d &edgemap);
    virtual EllipseMappingEdge2d* clone() const;
    virtual EllipseMappingEdge2d* inverse() const;
 
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    EllipseMappingEdge2d* part(const Real t0, const Real t1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Real2d m_;
    Real a_, b_;
    Real angle0_, angle1_;
 
    EllipseMappingEdge2d(const Real2d center, const Real a, const Real b,
                         const Real angle0, const Real angle1);
  };
 
  // ************************************************** ParabelMappingEdge2d **
 
  class ParabelMappingEdge2d : public MappingEdge2d {
  public:
    ParabelMappingEdge2d(const Real2d vtx0, const Real2d vtx1,
                         const Real2d cusp, const Real2d ap);
    ParabelMappingEdge2d(const ParabelMappingEdge2d &edgemap);
    virtual ParabelMappingEdge2d* clone() const;
 
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    ParabelMappingEdge2d* part(const Real t0, const Real t1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const UnitNd<2> v_, vrot_;
    const Real2d cusp_;
    const Real s0_, deltas_;
    Real p_;
    const Real s_(const Real t) const;
  };
 
  // ************************************************ MappingParallelEdge2d **
 
  class MappingParallelEdge2d : public MappingEdge2d {
  public:
    MappingParallelEdge2d(const MappingEdge2d &edgemap, const Real d);
    MappingParallelEdge2d(const MappingParallelEdge2d &edgemap);
    virtual MappingParallelEdge2d* clone() const;
 
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    virtual Real curvature(const Real t, const uint n = 0) const;
    virtual Real2d n0(const Real t) const;
 
    inline const MappingEdge2d& map() const { return *edgemap_; }
    inline const Real d() const { return d_; }
    MappingParallelEdge2d* part(const Real t0, const Real t1) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::unique_ptr<const MappingEdge2d> edgemap_;
    Real d_;
  };
 
  class MappingQuad2d;
 
  // ***************************************************** MappingQuadEdge2d **
 
  class MappingQuadEdge2d : public MappingEdge2d {
  public:
    MappingQuadEdge2d(const MappingQuad2d& map, const uint edge);
    MappingQuadEdge2d(const MappingQuad2d& map,
                      const Real2d xi0, const Real2d xi1);
    MappingQuadEdge2d(const MappingQuadEdge2d &edgemap);
    virtual MappingQuadEdge2d* clone() const;
 
    virtual Real2d operator()(const Real t) const;
    virtual Real2d derivative(const Real t, const uint n = 1) const;
    virtual MappingQuadEdge2d* inverse() const {
      return new MappingQuadEdge2d(*map_, xi_ + xidiff_, xi_);
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::unique_ptr<MappingQuad2d> map_;
    Real2d xi_, xidiff_;
 
    // Local coordinate of an point on the edge
    inline const Real2d lcoord_(const Real t) const {
      return xi_ + xidiff_*t;
    }
    // Derivative of the local mapping
    inline const Real2d lderiv_(const Real t) const {
      return xidiff_;
    }
  };
 
  // ***************************************************** MappingTriangle2d **
 
  class MappingTriangle2d : public Map2d {
  public:
    virtual Real2d operator()(const Real x, const Real y) const = 0;
 
    virtual MappingTriangle2d* clone() const = 0;
 
    //virtual bool curved() const = 0;
 
    virtual Real jacobianDeterminant(const Real x, const Real y) const;
 
    virtual MapReal2d jacobian(const Real x, const Real y) const = 0;
 
    virtual MapReal2d jacobianInverse(const Real x, const Real y) const;
  };
 
  // ****************************************************** VertexTriangle2d **
 
  class VertexTriangle2d : public MappingTriangle2d {
  public:
    VertexTriangle2d(Real2d vtx0, Real2d vtx1, Real2d vtx2);
 
    VertexTriangle2d(const VertexTriangle2d &map);
 
    virtual ~VertexTriangle2d();
    virtual VertexTriangle2d* clone() const;
    //virtual bool curved() const;
    
    virtual bool straight() const {return true;};
    
    virtual MapReal2d jacobian(const Real x, const Real y) const {
      return jacobian_;
    }
 
    Real2d operator()(const Real x, const Real y) const {
      conceptsAssert(0.0 <= y, Assertion());
      conceptsAssert(0.0 <= x, Assertion());
      conceptsAssert(y <= x, Assertion());
      conceptsAssert(x <= 1.0, Assertion());
 
      Real2d res(vtx1_);  res -= vtx0_;  res *= x;
      Real2d res2(vtx2_); res2 -= vtx1_; res2 *= y;
      res += res2;
 
      res += vtx0_;
      return res;
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Real2d vtx0_, vtx1_, vtx2_;
    MapReal2d jacobian_;
  };
 
  // ******************************************************* MapTriangle2d **
 
  class MapTriangle2d : public MappingTriangle2d {
  public:
    MapTriangle2d(char* map, const Real scX, const Real scY);
 
    MapTriangle2d(const MapTriangle2d &map);
 
    virtual ~MapTriangle2d();
    virtual MapTriangle2d* clone() const;
    //virtual bool curved() const;
    virtual MapReal2d jacobian(const Real x, const Real y) const;
 
    Real2d operator()(const Real x, const Real y) const {
      conceptsAssert(0.0 <= y, Assertion());
      conceptsAssert(0.0 <= x, Assertion());
      conceptsAssert(y <= x, Assertion());
      conceptsAssert(x <= 1.0, Assertion());
      return Real2d(map_, scx_ * x, scy_ * (x == 0 ? 0 : y/x));
    }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uchar* map_;
 
    uint sz_;
 
    const Real scx_;
 
    const Real scy_;
  };
 
 
  // ********************************************************* MappingQuad2d **
 
  class MappingQuad2d : public Map2d {
  public:
    virtual Real2d operator()(Real x, Real y) const = 0;
 
    virtual Real jacobianDeterminant(const Real x, const Real y) const;
 
    virtual MapReal2d jacobian(const Real x, const Real y) const = 0;
 
    virtual MapReal2d jacobianInverse(const Real x, const Real y) const;
 
    virtual MapReal2d hessian(uint i, const Real x, const Real y) const = 0;
 
    virtual MapReal2d inverseLaplace(const Real x, const Real y) const;
 
 
    virtual Real lineElement(const Real x, const uint edge) const;
 
    virtual MappingEdge2d* edge(const uint edge) const {
      return new MappingQuadEdge2d(*this, edge);
    }
 
    virtual MappingQuad2d* part(const Real2d x0, const Real2d y0) const;
 
    virtual MappingQuad2d* clone() const = 0;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
  // ************************************************************* MapQuad2d **
 
  class MapQuad2d : public MappingQuad2d {
  public:
    MapQuad2d(const char* map, Real scX, Real scY,
              Real2d org = Real2d(0.0, 0.0), Real stretch = 1.0);
    MapQuad2d(const char* map, Real scX, Real scY, Real2d org,
              Real2d scaling);
 
    MapQuad2d(const MapQuad2d &map);
    virtual ~MapQuad2d();
 
    virtual Real2d operator()(Real x, Real y) const { 
      Real2d v(map_, scx_ * x, scy_ * y);
      return Real2d(v[0]*scaling_[0] + org_[0], v[1]*scaling_[1] + org_[0]);
    }
 
    virtual MapReal2d jacobian(const Real x, const Real y) const;
 
    virtual MapReal2d hessian(uint i, const Real x, const Real y) const;
 
    virtual MapQuad2d* clone() const { return new MapQuad2d(*this); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    uchar* map_;
 
    uint sz_;
 
    Real scx_;
 
    Real scy_;
 
    Real2d org_;
    Real2d scaling_;
 
    void operator=(const MapQuad2d &);
  };
 
  // ********************************************************** VertexQuad2d **
 
  class VertexQuad2d : public MappingQuad2d {
  public:
    VertexQuad2d(Real2d vtx0, Real2d vtx1, Real2d vtx2, Real2d vtx3);
    VertexQuad2d(const VertexQuad2d& v);
    virtual ~VertexQuad2d();
    virtual Real2d operator()(Real x, Real y) const;
 
    virtual MapReal2d jacobian(const Real x, const Real y) const;
 
    virtual MapReal2d hessian(uint i, const Real x, const Real y) const;
 
    virtual Real lineElement(const Real x, const uint edge) const;
 
    virtual MappingStraightEdge2d* edge(const uint edge) const;
 
    virtual VertexQuad2d* part(const Real2d x0, const Real2d y0) const;
 
    virtual VertexQuad2d* clone() const;
 
    virtual bool straight() const {return true;};
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Real2d vtx_[4];
    mutable Real len_[4];
    bool hasLength_;
  };
 
  // ********************************************************* BlendingQuad2d **
 
  class BlendingQuad2d : public MappingQuad2d {
  public:
    BlendingQuad2d(const MappingEdge2d* edgemap0,
                   const MappingEdge2d* edgemap1,
                   const MappingEdge2d* edgemap2,
                   const MappingEdge2d* edgemap3 = 0,
                   const Real2d* vtx0 = 0, const Real2d* vtx1 = 0,
                   const Real2d* vtx2 = 0, const Real2d* vtx3 = 0);
    BlendingQuad2d(const BlendingQuad2d& v);
    ~BlendingQuad2d();
    virtual Real2d operator()(Real x, Real y) const;
    virtual MapReal2d jacobian(const Real x, const Real y) const;
    virtual MapReal2d hessian(uint i, const Real x, const Real y) const;
    virtual Real lineElement(const Real x, const uint edge) const;
    virtual MappingEdge2d* edge(const uint edge) const {
      conceptsAssert(edge < 4, Assertion());
      return edgemap_[edge]->clone();
    }
    virtual BlendingQuad2d* clone() const { return new BlendingQuad2d(*this); }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Real2d vtx_[4];
    const MappingEdge2d* edgemap_[4];
    void testJacobian_() const;
  };
 
  // ***************************************************** PartMappingQuad2d **
 
  class PartMappingQuad2d : public MappingQuad2d {
  public:
    PartMappingQuad2d(const MappingQuad2d& map, 
                      const Real2d xi0, const Real2d xi1);
    PartMappingQuad2d(const PartMappingQuad2d &map);
    virtual PartMappingQuad2d* clone() const;
    virtual Real2d operator()(Real x, Real y) const;
    virtual MapReal2d jacobian(const Real x, const Real y) const;
    virtual MapReal2d hessian(uint i, const Real x, const Real y) const;
    virtual MappingEdge2d* edge(const uint edge) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    std::unique_ptr<const MappingQuad2d> map_;
    const Real2d x0_;
    const Real2d d_;
 
    inline const Real2d xi_(const Real x, const Real y) const {
      return x0_ + Real2d(x * d_[0], y * d_[1]);
    }
  };
 
  // ********************************************* InverseVertexQuadSector2d **
 
  class InverseVertexQuadSector2d : public Map2d {
  public: 
    InverseVertexQuadSector2d(const Real2d vtx0, const Real2d vtx1,
                              const Real2d vtx2, const Real2d vtx3);
    InverseVertexQuadSector2d(const Sequence<Real2d> x);
    // Copy constructor
    InverseVertexQuadSector2d(const InverseVertexQuadSector2d& v);
    virtual ~InverseVertexQuadSector2d();
    virtual Real2d operator()(const Real2d x) const;
    virtual MapReal2d jacobian(const Real2d x) const;
    virtual bool is_inregion(const Real2d x) const;
    virtual bool is_inpatch(const Real2d x) const;
    virtual InverseVertexQuadSector2d* clone() const;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    
  private:
    void construct_(const Real2d vtx0, const Real2d vtx1, 
                    const Real2d vtx2, const Real2d vtx3);
    Real2d sector_p_[4];
    Real2d vtx_[4]; 
  };
 
  // ********************************************************* MapTriangle3d **
 
  class MapTriangle3d : public Map2d {
  public:
    MapTriangle3d(const char* map, Real scX, Real scY);
 
    MapTriangle3d(Real3d vtx0, Real3d vtx1, Real3d vtx2);
 
    MapTriangle3d(const MapTriangle3d &map);
    ~MapTriangle3d() { delete[] map_; };
 
    Real3d operator()(Real x, Real y) const { 
      if (map_)
        return Real3d(map_, scx_ * x, scy_ * y);
      else {
        Real3d res;
        res.lincomb(B1_, B2_, x, x*y);
        res += b_;
        return res;
      }
    }
 
    inline MapTriangle3d* clone() const {return new MapTriangle3d(*this);};
  private:
    uchar* map_;
 
    uint sz_;
 
    Real scx_;
 
    Real scy_;
 
    Real3d B1_, B2_, b_;
 
    void operator=(const MapTriangle3d &);
  };
 
} // namespace concepts
 
#endif // elementMaps_hh