arrays.hh

Go to the documentation of this file.

 
#include "toolbox/array.hh"
#include "cell1D.hh"
#include "cell2D.hh"
#include "cell3D.hh"
 
#ifndef geometryArrays_hh
#define geometryArrays_hh
 
namespace concepts {
 
  // ************************************************************* makeArray **
 
  template<class F, class G>
  void makeArray(const F& cell, const Array<Real>& p, G (F::*fun)(Real) const,
      Array<G>& array) {
    array.resize(p.size());
    G* coord = array;
    const Real* pAcc = p;
    for (uint i = 0; i < p.size(); ++i)
      *coord++ = (cell.*fun)(*pAcc++);
  }
 
  template<class F, class G>
  void makeArray(const F& cell, const Array<Real>& pX, const Array<Real>& pY,
                 G (F::*fun)(Real, Real) const, Array<G>& array,
                 bool istensor = true)
  {
    if (istensor) {
      array.resize(pX.size() * pY.size());
      G* coord = array;
      const Real* pxAcc = pX;
      for (uint i = 0; i < pX.size(); ++i) {
        const Real* pyAcc = pY;
        for (uint j = 0; j < pY.size(); ++j)
          *coord++ = (cell.*fun)(*pxAcc, *pyAcc++);
        pxAcc++;
      }
    } else { // non tensor case, simple one loop over all points
      conceptsAssert(pX.size() == pY.size(), concepts::Assertion());
 
      array.resize(pX.size());
      G* coord = array;
      const Real* pxAcc = pX;
      const Real* pyAcc = pY;
      for (uint i = 0; i < pX.size(); ++i)
        *coord++ = (cell.*fun)(*pxAcc++, *pyAcc++);
    } // else
  }
 
  template<class F, class G>
  void makeArray(const F& cell, 
                 const Array<Real>& pX,
                 const Array<Real>& pY,
                 const Array<Real>& pZ,
                 G (F::*fun)(Real, Real, Real) const, Array<G>& array,
                 bool istensor = true)
  {
    if (istensor) {
      array.resize(pX.size() * pY.size() * pZ.size());
      G* coord = array;
      const Real* pxAcc = pX;
      for (uint i = 0; i < pX.size(); ++i) {
        const Real* pyAcc = pY;
        for (uint j = 0; j < pY.size(); ++j) {
          const Real* pzAcc = pZ;
          for (uint k = 0; k < pZ.size(); ++k)
            *coord++ = (cell.*fun)(*pxAcc, *pyAcc, *pzAcc++);
          pyAcc++;
        }
        pxAcc++;
      }
    } else { // non tensor case, simple one loop over all points
      conceptsAssert(pX.size() == pY.size(), concepts::Assertion());
      conceptsAssert(pX.size() == pZ.size(), concepts::Assertion());
 
      array.resize(pX.size());
      G* coord = array;
      const Real* pxAcc = pX;
      const Real* pyAcc = pY;
      const Real* pzAcc = pZ;
      for (uint i = 0; i < pX.size(); ++i)
        *coord++ = (cell.*fun)(*pxAcc++, *pyAcc++, *pzAcc++);
    } // else
  }
 
  // ************************************************************ ArrayCoord **
 
  template<uint dim>
  class ArrayCoord {
  };
 
  // ********************************************************* ArrayCoord<1> **
 
  template<>
  class ArrayCoord<1> : public Array<Real> {
  public:
    ArrayCoord(const Edge1d& edge, const Array<Real>& p);
  };
 
  // ********************************************************* ArrayCoord<2> **
 
  template<>
  class ArrayCoord<2> : public Array<Real2d> {
  public:
    ArrayCoord(const Quad2d& quad,
               const Array<Real>& qX,
               const Array<Real>& qY,
               bool tensor = true);
    ArrayCoord(const Edge2d& edge, const Array<Real>& p);
  };
 
  // ********************************************************* ArrayJacobian **
 
  template<int gdim, int ldim>
  class ArrayJacobian {
  };
 
  // **************************************************** ArrayJacobian<1,1> **
 
  template<>
  class ArrayJacobian<1, 1> : public Array<Real> {
  public:
    ArrayJacobian(const EdgeNd& edge, const Array<Real>& p);
    ArrayJacobian(const QuadNd& quad,
                  const Array<Real>& qX,
                  const Array<Real>& qY,
                  bool tensor = true);
  };
 
  // **************************************************** ArrayJacobian<2,2> **
 
  template<>
  class ArrayJacobian<2, 2> : public Array<Mapping<Real, 2> > {
  public:
    ArrayJacobian(const QuadNd& quad,
                  const Array<Real>& qX,
                  const Array<Real>& qY,
                  bool istensor = true);
 
    ArrayJacobian(const QuadNd& quad, uint k, const Array<Real>& q);
  };
 
  // **************************************************** ArrayJacobian<2,2> **
 
  template<>
  class ArrayJacobian<3, 3> : public Array<Mapping<Real, 3> > {
  public:
    ArrayJacobian(const Hexahedron3d& quad,
                  const Array<Real>& qX,
                  const Array<Real>& qY,
                  const Array<Real>& qZ,
                  bool istensor = true);
  };
 
  // ******************************************************** ArrayLocalCoord **
 
  class ArrayLocalCoord: public Array<Real> {
  public:
    ArrayLocalCoord(const QuadNd& quad, uint k, uint dim, const Array<Real>& q);
  private:
    inline bool alongEdge(uint k, uint dim) {
      return (k == 0 && dim == 1) || (k == 3 && dim == 0)
          || (k == 1 && dim == 0) || (k == 2 && dim == 1);
    }
  };
 
 
   template<int gdim, int ldim>
   class ArrayHessian {
   };
 
   // **************************************************** ArrayHessian<1,1> **
 
   template<>
   class ArrayHessian<1, 1> : public Array<Real> {
   public:
     ArrayHessian(const EdgeNd& edge, const Array<Real>& p);
 
   };
 
 
} // namespace concepts
 
#endif // geometryArrays_hh