geometry/chart.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_CHART_H
#define _MATHTL_CHART_H

#include <iostream>
#include <string>
#include <algebra/matrix.h>
#include <algebra/vector.h>
#include <geometry/point.h>
#include <utils/fixed_array1d.h>

using std::string;

namespace MathTL
{
  template <unsigned int DIM_d, unsigned int DIM_m = DIM_d>
  class Chart
  {
  public:
    virtual ~Chart() {}

    virtual void map_point(const Point<DIM_d>& x, Point<DIM_m>& y) const = 0;

    virtual double map_point(const double, const int) const = 0;

    virtual void map_point_inv(const Point<DIM_m>& x, Point<DIM_d>& y) const = 0;

    virtual double map_point_inv(const double, const int) const = 0;

    virtual const double Gram_factor(const Point<DIM_d>& x) const = 0;

    virtual const double Gram_D_factor(const unsigned int i,
                                       const Point<DIM_d>& x) const = 0;

    virtual const double Dkappa_inv(const unsigned int i,
                                    const unsigned int j,
                                    const Point<DIM_d>& x) const = 0;

    virtual const double a_i(const int i) const = 0;
    
    virtual const bool in_patch(const Point<DIM_m>& x) const = 0;

    virtual const string to_string() const = 0;
  };
  
  template <unsigned int DIM_d, unsigned int DIM_m>
  std::ostream& operator << (std::ostream& s, const Chart<DIM_d,DIM_m>&);

  //
  // Some examples:

  template <unsigned int DIM>
  class AffineLinearMapping
    : public Chart<DIM,DIM>
  {
  public:
    AffineLinearMapping();

    AffineLinearMapping(const Matrix<double>& A, const Point<DIM>& b);

    virtual ~AffineLinearMapping() {}

    void map_point(const Point<DIM>&, Point<DIM>&) const;
    double map_point(const double, const int) const;

    void map_point_inv(const Point<DIM>&, Point<DIM>&) const;
    double map_point_inv(const double, const int) const;

    const double Gram_factor(const Point<DIM>&) const;
    const double Gram_D_factor(const unsigned int i, const Point<DIM>& x) const;
    const double Dkappa_inv(const unsigned int i, const unsigned int j,
                            const Point<DIM>& x) const;
    const bool in_patch(const Point<DIM>& x) const;
    
    const double a_i(const int i) const { return A_(i,i); };


    const string to_string() const;

    static const string className;

    const Matrix<double>& A() const { return A_; }
    
    const Point<DIM>& b() const { return b_; }

  protected:
    Matrix<double> A_, A_inv;
    double det_A;
    double square_root_of_abs_det_A;
    Point<DIM> b_;
  };

  template <unsigned int DIM>
  class SimpleAffineLinearMapping
    : public Chart<DIM,DIM>
  {
    public:
    virtual ~SimpleAffineLinearMapping() {};


    SimpleAffineLinearMapping();

    SimpleAffineLinearMapping(const FixedArray1D<double,DIM>& A, const Point<DIM>& b);

    void map_point(const Point<DIM>&, Point<DIM>&) const;
    double map_point(const double, const int) const {return 0.;};// dummy;

    void map_point_inv(const Point<DIM>&, Point<DIM>&) const;
    double map_point_inv(const double, const int) const {return 0.;};// dummy

    const double Gram_factor(const Point<DIM>&) const;
    const double Gram_D_factor(const unsigned int i, const Point<DIM>& x) const;
    const double Dkappa_inv(const unsigned int i, const unsigned int j,
                            const Point<DIM>& x) const;
    const bool in_patch(const Point<DIM>& x) const;

    const double a_i(const int i) const { return A_[i]; };

    const string to_string() const;

    static const string className;

    const FixedArray1D<double,DIM>& A() const { return A_; }
    
    const Point<DIM>& b() const { return b_; }

  protected:
    FixedArray1D<double,DIM> A_, A_inv;
    double det_A;
    double square_root_of_abs_det_A;
    Point<DIM> b_;
    
    
  };

  class LinearBezierMapping : public Chart<2,2>
  {

  public:
    LinearBezierMapping ();

    LinearBezierMapping (const LinearBezierMapping&);

    LinearBezierMapping (const Point<2> &, const Point<2> &,
                         const Point<2> &, const Point<2> &);

    ~LinearBezierMapping () {};
    
    LinearBezierMapping& operator = (const LinearBezierMapping& x);

    void setup ();

    const Point<2>& get_b_00() const;
    const Point<2>& get_b_01() const;
    const Point<2>& get_b_10() const;
    const Point<2>& get_b_11() const;
 
    void map_point(const Point<2>&, Point<2>&) const;
    double map_point(const double, const int) const {return 0.;};// dummy;

    void map_point_inv(const Point<2>&, Point<2>&) const;
    double map_point_inv(const double, const int) const {return 0.;};// dummy;

    const double Gram_factor(const Point<2>& x) const;
    const double Gram_D_factor(const unsigned int i, const Point<2>& x) const;
    const double Dkappa_inv(const unsigned int i, const unsigned int j,
                            const Point<2>& x) const;
    const bool in_patch(const Point<2>& x) const;

    const double a_i(const int i) const { return 0.; };
    
    const string to_string() const;
   
    static const string className;

    //vertices of the quadrangle to parametrize
    Point<2> b_00;    
    Point<2> b_01;
    Point<2> b_10;
    Point<2> b_11;

  private:
    //the following two routines are helpers that are only of interest
    //in this special 2D BezierMapping case!
    //only for the special case of LinearBezierMapping
    const double det_DKappa(const Point<2>& p) const;

    //partial derivative with respect to i-th component
    //of j-th component of Kappa
    const double partial_i_Kappa_j(const unsigned int i,
                                   const unsigned int j,
                                   const Point<2>& x) const;

    //vertices of generic qudrangle, needed for inverting this mapping,
    //initialized in constructor
    Point<2> b_gen_00;
    Point<2> b_gen_10;
    Point<2> b_gen_01;
    Point<2> b_gen_11;

    Vector<double> d_ds_d_dt_kappa_r;
    Vector<double> d_dt_d_ds_kappa_r;

    Vector<double> min_b00_plus_b10;
    Vector<double> min_b00_plus_b01;

    //some quantities, worth storing to prevent dispensable recomutation
    double cos_rot_angle;
    double sin_rot_angle;
    double rot_angle;
    double shearing_param;
    double scaleX;
    double scaleY;
    
    //signum of det( D (mapPoint(x,y) ) )
    //identical for all (x,y)!
    bool sgn_det_D;


  };
}

#include "geometry/chart.cpp"

#endif