geometry/elementMaps.hh

Go to the documentation of this file.

/* Element maps.
 * Geometrical information (coordinates)
 */

#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 << "Map1d()";};
  };

  // ***************************************************************** Map2d **

  class Map2d : public OutputOperator {
  public:
    std::ostream& info(std::ostream& os) const {return os << "Map2d()";};
    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::auto_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;

    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::auto_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::auto_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 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);

    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];
    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::auto_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