elementMaps3D.hh

Go to the documentation of this file.

 
#ifndef elementMaps3D_hh
#define elementMaps3D_hh
 
#include "basics/outputOperator.hh"
#include "basics/typedefs.hh"
#include "basics/vectorsMatrices.hh"
#include "basics/debug.hh"
#include "elementMaps.hh"
 
// debugging
#define geoMapTetrahedron3dAppl_D 0
 
namespace concepts {
 
  // ***************************************************************** Map3d **
 
  class Map3d : public OutputOperator {
  public:
    std::ostream& info(std::ostream& os) const {return os << "Map3d()";};
  };
 
 
  // ****************************************************** MapTetrahedron3d **
 
  class MapTetrahedron3d : public Map3d {
  public:
    MapTetrahedron3d(char* map, Real scX, Real scY, Real scZ);
 
    MapTetrahedron3d(Real3d vtx0, Real3d vtx1, Real3d vtx2, Real3d vtx3);
 
    MapTetrahedron3d(const MapTetrahedron3d &map);
    ~MapTetrahedron3d() { delete[] map_; }
 
    Real3d operator()(Real x, Real y, Real z) const { 
      if (map_)
  return Real3d(map_, scx_ * x, scy_ * y, scz_ * z);
      else {
  Real3d xi(x*(1.0-y)*(1.0-z), y*(1.0-z), z);
  Real3d res(B_ * xi);
  res += b_;
  DEBUGL(geoMapTetrahedron3dAppl_D, '(' << x << ", " << y 
         << ", " << z << ") -> " << xi << " -> " << res);
  return res;
      }
    }
 
    Real jacobian() const {
      if (map_)
        throw conceptsException(MissingFeature("jacobian not supported"));
      return B_.determinant(); }
 
    inline MapTetrahedron3d* clone() const {
      return new MapTetrahedron3d(*this);
    };
 
private:
    uchar* map_;
 
    uint sz_;
 
    Real scx_;
 
    Real scy_;
 
    Real scz_;
 
    MapReal3d B_;
 
    Real3d b_;
 
    void operator=(const MapTetrahedron3d &);
  };
  
  
    // ********************************************************** MappingHexahedron3d **
 
  //forward declaration
  class MappingQuad3d;  
  
  class MappingHexahedron3d : public Map3d  {
  public:
    
    virtual ~MappingHexahedron3d() {}
    virtual Real3d operator()(Real x, Real y, Real z) const = 0;
    
    virtual MapReal3d jacobian(const Real x, const Real y, const Real z) const = 0;
    
    virtual MapReal3d hessian(uint i,const Real x, const Real y, const Real z) const = 0;
    
    virtual MapReal3d jacobianInverse(const Real x, const Real y, const Real z) const {
        return jacobian(x,y,z).inverse(); }
        
    virtual Real jacobianDeterminant(const Real x, const Real y, const Real z) const {
        return jacobian(x,y,z).determinant(); }
        
    virtual MappingQuad3d* face(const uint face) const = 0;
    
    virtual MappingHexahedron3d* clone() const = 0;
 
    virtual MappingHexahedron3d* part(const Real3d xi0, const Real3d xi1) const = 0;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const = 0;
  };
 
  // ********************************************************* BlendingHexahedron3d **
 
  //forward declaration
  class MappingQuad3d;
  class MappingEdge3d;
  class Hexahedron;
  
  class BlendingHexahedron3d : public MappingHexahedron3d {
  public:
 
    BlendingHexahedron3d(const Sequence<const MappingQuad3d*> quadmap, const Hexahedron& hex);
    BlendingHexahedron3d(const MappingQuad3d* quadmap0, const MappingQuad3d* quadmap1,
                         const MappingQuad3d* quadmap2, const MappingQuad3d* quadmap3,
                         const MappingQuad3d* quadmap4, const MappingQuad3d* quadmap5,
                         const Hexahedron& hex);
    BlendingHexahedron3d(const BlendingHexahedron3d& v);
    BlendingHexahedron3d* clone() const override { return new BlendingHexahedron3d(*this); }
    ~BlendingHexahedron3d() override;
    Real3d operator()(Real x, Real y, Real z) const override;
    MapReal3d jacobian(const Real x, const Real y, const Real z) const override;
    MapReal3d hessian(uint i,const Real x, const Real y, const Real z) const override;
    MappingQuad3d* face(const uint face) const override;
    MappingEdge3d* edge(const uint edge) const;
    MappingHexahedron3d* part(const Real3d xi0, const Real3d xi1) const override;
  protected:
    std::ostream& info(std::ostream& os) const override;
  private:
    // Maps of the faces
    std::array<std::unique_ptr<MappingQuad3d>,6> quadmap_;
    // Maps of the edges
    std::array<std::unique_ptr<MappingEdge3d>,12> edgemap_;
    // Coordinates of the vertices
    std::array<Real3d,8> vtx_;
 
    // Find all geometric boundary quads, edges and vertices
    // This function should actually be a member of a MappingHex class, which should independent of the
    // dimension of the ambient space    
    void construct_(const std::array<const MappingQuad3d*,6>& quadmap, const Hexahedron& hex);
  };
 
  // ***************************************************** PartMappingHexahedron3d **
  
  class PartMappingHexahedron3d : public MappingHexahedron3d {
  public:
    PartMappingHexahedron3d(const MappingHexahedron3d& map, 
                            const Real3d x0, const Real3d x1);
    PartMappingHexahedron3d(const PartMappingHexahedron3d& other);
    ~PartMappingHexahedron3d() override {};
    PartMappingHexahedron3d* clone() const override;
    Real3d operator()(Real x, Real y, Real z) const override;
    MapReal3d jacobian(const Real x, const Real y, const Real z) const override;
    MapReal3d hessian(uint i,const Real x, const Real y, const Real z) const override;
    MappingQuad3d* face(const uint face) const override;
    PartMappingHexahedron3d* part(const Real3d x0, const Real3d x1) const override;    
  protected:
    std::ostream& info(std::ostream& os) const override;
  private:
    std::unique_ptr<MappingHexahedron3d> map_;
    const Real3d x0_;
    const Real3d d_;
 
    inline const Real3d xi_(const Real x, const Real y, const Real z) const {
      return x0_ + Real3d(x * d_[0], y * d_[1], z * d_[2]);
    }
  };
 
 
  // ******************************************************* MapHexahedron3d **
 
  class MapHexahedron3d : public MappingHexahedron3d {
 
    uchar* map_;
 
    uint sz_;
 
    Real scx_;
 
    Real scy_;
 
    Real scz_;
 
    Real3d *comp_;
 
    void operator=(const MapHexahedron3d &);
 
  public:
    MapHexahedron3d(char* map, Real scX, Real scY, Real scZ);
 
    MapHexahedron3d(const Real3d& vtx0, const Real3d& vtx1, 
        const Real3d& vtx2, const Real3d& vtx3,
        const Real3d& vtx4, const Real3d& vtx5,
        const Real3d& vtx6, const Real3d& vtx7);
 
    MapHexahedron3d(const MapHexahedron3d &map);
    ~MapHexahedron3d() { delete[] map_;  delete[] comp_; }
 
    inline Real3d operator()(Real x, Real y, Real z) const { 
      conceptsAssert(x >= 0.0, Assertion());
      conceptsAssert(x <= 1.0, Assertion());
      conceptsAssert(y >= 0.0, Assertion());
      conceptsAssert(y <= 1.0, Assertion());
      conceptsAssert(z >= 0.0, Assertion());
      conceptsAssert(z <= 1.0, Assertion());
      if (map_)
  return Real3d(map_, scx_ * x, scy_ * y, scz_* z);
      else {
  Real3d res(comp_[0]*(x*y*z));
  res += comp_[1]*(x*y);
  res += comp_[2]*(x*z);
  res += comp_[3]*(y*z);
  res += comp_[4]*x;
  res += comp_[5]*y;
  res += comp_[6]*z;
  res += comp_[7];
  return res;
      }
    }
 
    MapReal3d jacobian(const Real x, const Real y, const Real z) const;
    
    MapReal3d hessian(uint i,const Real x, const Real y, const Real z) const;
 
    MapReal3d jacobianInverse(const Real x, const Real y, const Real z) const {
      return jacobian(x, y, z).inverse();
    }
 
    Real jacobianDeterminant(const Real x, const Real y, const Real z) const {
      return jacobian(x, y, z).determinant();
    }
 
    MapHexahedron3d* clone() const { return new MapHexahedron3d(*this); }
    
    MappingQuad3d* face(uint face) const;
    
    MapHexahedron3d* part(const Real3d x0, const Real3d z0) const;
    
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  };
 
 
  // *************************************************** MapParallelepiped3d **
 
  class MapParallelepiped3d : public Map3d {
  public:
    MapParallelepiped3d(Real3d vtx0, Real3d vtx1, Real3d vtx2, Real3d vtx3);
 
    MapParallelepiped3d(const MapParallelepiped3d &map);
 
    Real3d operator()(Real x, Real y, Real z) const {
      Real3d res(B1_*x + B2_*y + B3_*z);
      res += b_;
      return res;
    }
 
    Real jacobian() const {
      return (B1_ ^ B2_) * B3_;
    }
 
    inline MapParallelepiped3d* clone() const {
      return new MapParallelepiped3d(*this);
    };
 
  private:
    Real3d b_;
    Real3d B1_;
    Real3d B2_;
    Real3d B3_;
 
    void operator=(const MapParallelepiped3d &);
  };
 
} // namespace concepts
 
#endif // elementMaps3D_hh