[cig-commits] commit: Add a skeletal Tensor2_antisymmetric
Mercurial
hg at geodynamics.org
Fri Mar 9 13:25:00 PST 2012
changeset: 8:af762f2efe25
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Fri Mar 09 11:11:37 2012 -0800
files: FTensor.hpp Tensor2_antisymmetric.hpp Tensor2_antisymmetric/Tensor2_antisymmetric_Expr.hpp Tensor2_antisymmetric/Tensor2_antisymmetric_constructor.hpp Tensor2_antisymmetric/Tensor2_antisymmetric_value.hpp
description:
Add a skeletal Tensor2_antisymmetric
diff -r 8c683302bc11 -r af762f2efe25 FTensor.hpp
--- a/FTensor.hpp Fri Mar 09 09:32:56 2012 -0800
+++ b/FTensor.hpp Fri Mar 09 11:11:37 2012 -0800
@@ -27,6 +27,10 @@ namespace FTensor {
template <class T, int Dim> class Tensor2_symmetric;
template<class A, class T, int Dim, char i, char j>
class Tensor2_symmetric_Expr;
+
+ template <class T, int Dim> class Tensor2_antisymmetric;
+ template<class A, class T, int Dim, char i, char j>
+ class Tensor2_antisymmetric_Expr;
template <class A, class T, int Dim0, int Dim1, int Dim2,
char i, char j, char k> class Tensor3_Expr;
@@ -69,6 +73,7 @@ namespace FTensor {
#include "Tensor1.hpp"
#include "Tensor2.hpp"
#include "Tensor2_symmetric.hpp"
+#include "Tensor2_antisymmetric.hpp"
#include "Tensor3/Tensor3_Expr.hpp"
#include "Tensor3/Tensor3_contracted.hpp"
#include "Tensor3_dg.hpp"
diff -r 8c683302bc11 -r af762f2efe25 Tensor2_antisymmetric.hpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Tensor2_antisymmetric.hpp Fri Mar 09 11:11:37 2012 -0800
@@ -0,0 +1,10 @@
+/* Declarations for a Tensor2_antisymmetric. This only has a single
+ dimension parameter, since the dimensions of the two indices must
+ be the same for it to be antisymmetric. */
+
+#include "Tensor2_antisymmetric/Tensor2_antisymmetric_constructor.hpp"
+#include "Tensor2_antisymmetric/Tensor2_antisymmetric_value.hpp"
+// TODO implement Tensor2_antisymmetric_pointer.hpp
+
+
+#include "Tensor2_antisymmetric/Tensor2_antisymmetric_Expr.hpp"
diff -r 8c683302bc11 -r af762f2efe25 Tensor2_antisymmetric/Tensor2_antisymmetric_Expr.hpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Tensor2_antisymmetric/Tensor2_antisymmetric_Expr.hpp Fri Mar 09 11:11:37 2012 -0800
@@ -0,0 +1,80 @@
+/* Declares a wrapper class for antisymmetric rank 2 Tensor expressions.
+ It is specialized for Tensor3_number_rhs_2. */
+
+template<class A, class T, int Dim, char i, char j>
+class Tensor2_antisymmetric_Expr
+{
+ A iter;
+public:
+ Tensor2_antisymmetric_Expr(A &a): iter(a) {}
+ T operator()(const int N1, const int N2) const
+ {
+ return iter(N1,N2);
+ }
+};
+
+template<class A, class T, int Tensor_Dim, int Dim, char i, char j>
+class Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j>
+{
+ Tensor2_antisymmetric<A,Tensor_Dim> &iter;
+public:
+ Tensor2_antisymmetric_Expr(Tensor2_antisymmetric<A,Tensor_Dim> &a): iter(a) {}
+ T & operator()(const int N1, const int N2)
+ {
+ return iter(N1,N2);
+ }
+ T operator()(const int N1, const int N2) const
+ {
+ return iter(N1,N2);
+ }
+
+ // /* Various assignment operators. I have to explicitly declare the
+ // second operator= because otherwise the compiler will generate its
+ // own and not use the template code. */
+
+ // template<class B, class U>
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator=(const Tensor2_antisymmetric_Expr<B,U,Dim,i,j> &result);
+
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator=(const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &result);
+
+ // template<class B, class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator+=(const Tensor2_antisymmetric_Expr<B,U,Dim,i,j> &result);
+
+ // template<class B, class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator-=(const Tensor2_antisymmetric_Expr<B,U,Dim,i,j> &result);
+
+ // template<class B, class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator&=(const Tensor2_antisymmetric_Expr<B,U,Dim,i,j> &result);
+
+ // /* This is for when the indices are switched (i,j) -> (j,i). */
+
+ // template<class B, class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator=(const Tensor2_antisymmetric_Expr<B,U,Dim,j,i> &result);
+
+ // template<class B, class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator+=(const Tensor2_antisymmetric_Expr<B,U,Dim,j,i> &result);
+
+ // template<class B, class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> &
+ // operator-=(const Tensor2_antisymmetric_Expr<B,U,Dim,j,i> &result);
+
+ // /* Operations with just generics. */
+
+ // template<class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> & operator=(const U &d);
+ // template<class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> & operator+=(const U &d);
+ // template<class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> & operator-=(const U &d);
+ // template<class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> & operator*=(const U &d);
+ // template<class U> inline
+ // const Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<A,Tensor_Dim>,T,Dim,i,j> & operator/=(const U &d);
+};
diff -r 8c683302bc11 -r af762f2efe25 Tensor2_antisymmetric/Tensor2_antisymmetric_constructor.hpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Tensor2_antisymmetric/Tensor2_antisymmetric_constructor.hpp Fri Mar 09 11:11:37 2012 -0800
@@ -0,0 +1,43 @@
+/* A helper class that allows simple initialization of the Tensor2_antisymmetric,
+ but only if it has the correct number of elements. */
+
+template<class T, int Tensor_Dim>
+class Tensor2_antisymmetric_constructor;
+
+template<class T>
+class Tensor2_antisymmetric_constructor<T,2>
+{
+public:
+ Tensor2_antisymmetric_constructor(T data[(1*2)/2], T d01)
+ {
+ data[0]=d01;
+ }
+};
+
+template<class T>
+class Tensor2_antisymmetric_constructor<T,3>
+{
+public:
+ Tensor2_antisymmetric_constructor(T data[(2*3)/2], T d01, T d02, T d12)
+ {
+ data[0]=d01;
+ data[1]=d02;
+ data[2]=d12;
+ }
+};
+
+template<class T>
+class Tensor2_antisymmetric_constructor<T,4>
+{
+public:
+ Tensor2_antisymmetric_constructor(T data[(3*4)/2], T d01, T d02, T d03,
+ T d12, T d13, T d23)
+ {
+ data[0]=d01;
+ data[1]=d02;
+ data[2]=d03;
+ data[3]=d12;
+ data[4]=d13;
+ data[5]=d23;
+ }
+};
diff -r 8c683302bc11 -r af762f2efe25 Tensor2_antisymmetric/Tensor2_antisymmetric_value.hpp
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Tensor2_antisymmetric/Tensor2_antisymmetric_value.hpp Fri Mar 09 11:11:37 2012 -0800
@@ -0,0 +1,197 @@
+/* A general version, not for pointers. */
+
+template <class T, int Tensor_Dim>
+class Tensor2_antisymmetric
+{
+ T data[(Tensor_Dim*(Tensor_Dim-1))/2];
+public:
+ Tensor2_antisymmetric() {}
+
+ /* Tensor_Dim=2 */
+ Tensor2_antisymmetric(T d01)
+ {
+ Tensor2_antisymmetric_constructor<T,Tensor_Dim>(data,d01);
+ }
+
+ /* Tensor_Dim=3 */
+ Tensor2_antisymmetric(T d01, T d02, T d12)
+ {
+ Tensor2_antisymmetric_constructor<T,Tensor_Dim>(data,d01,d02,d12);
+ }
+
+ /* Tensor_Dim=4 */
+ Tensor2_antisymmetric(T d01, T d02, T d03, T d12, T d13, T d23)
+ {
+ Tensor2_antisymmetric_constructor<T,Tensor_Dim>
+ (data,d01,d02,d03,d12,d13,d23);
+ }
+
+ /* There are two ways of accessing the values inside,
+ unsafe(int,int) and operator(int,int). unsafe(int,int) will give
+ you a wrong answer if you aren't careful. The problem is that we
+ only store the minimal set of components, but some have different
+ signs. We can't return the negative of a component, and assign
+ something to it, because that would assign something to a
+ temporary. To get the correct answer if you don't want to change
+ the value, just use operator(int,int). */
+
+ T & unsafe(const int N1, const int N2)
+ {
+#ifdef FTENSOR_DEBUG
+ if(N1>=Tensor_Dim || N1<0 || N2>=Tensor_Dim || N2<0 || N1>=N2)
+ {
+ std::stringstream s;
+ s << "Bad index in Tensor2_antisymmetric<T," << Tensor_Dim
+ << ">.operator(" << N1 << "," << N2 << ")"
+ << std::endl;
+ throw std::runtime_error(s.str());
+ }
+#endif
+ return data[(N2-1)+(N1*(2*(Tensor_Dim-1)-N1-1))/2];
+ }
+
+ T operator()(const int N1, const int N2) const
+ {
+#ifdef FTENSOR_DEBUG
+ if(N1>=Tensor_Dim || N1<0 || N2>=Tensor_Dim || N2<0)
+ {
+ std::stringstream s;
+ s << "Bad index in Tensor2_antisymmetric<T," << Tensor_Dim
+ << ">.operator(" << N1 << "," << N2 << ") const"
+ << std::endl;
+ throw std::runtime_error(s.str());
+ }
+#endif
+ return N1==N2 ? 0 : (N1<N2 ? data[(N2-1)+(N1*(2*(Tensor_Dim-1)-N1-1))/2]
+ : -data[(N1-1)+(N2*(2*(Tensor_Dim-1)-N2-1))/2]);
+ }
+
+ /* These operator()'s are the first part in constructing template
+ expressions. They can be used to slice off lower dimensional
+ parts. They are not entirely safe, since you can accidently use a
+ higher dimension than what is really allowed (like Dim=5). */
+
+ /* This returns a Tensor2_Expr, since the indices are not really
+ antisymmetric anymore since they cover different dimensions. */
+
+ template<char i, char j, int Dim0, int Dim1>
+ Tensor2_Expr<Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j>
+ operator()(const Index<i,Dim0> index1, const Index<j,Dim1> index2)
+ {
+ return Tensor2_Expr<Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j>
+ (*this);
+ }
+
+ template<char i, char j, int Dim0, int Dim1>
+ Tensor2_Expr<const Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j>
+ operator()(const Index<i,Dim0> index1, const Index<j,Dim1> index2) const
+ {
+ return Tensor2_Expr<const Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim0,Dim1,i,j>
+ (*this);
+ }
+
+ /* This returns a Tensor2_antisymmetric_Expr, since the indices are still
+ antisymmetric on the lower dimensions. */
+
+ template<char i, char j, int Dim>
+ Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim,i,j>
+ operator()(const Index<i,Dim> index1, const Index<j,Dim> index2)
+ {
+ return Tensor2_antisymmetric_Expr<Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim,i,j>
+ (*this);
+ }
+
+ template<char i, char j, int Dim>
+ Tensor2_antisymmetric_Expr<const Tensor2_antisymmetric<T,Tensor_Dim>,T,Dim,i,j>
+ operator()(const Index<i,Dim> index1, const Index<j,Dim> index2) const
+ {
+ return Tensor2_antisymmetric_Expr<const Tensor2_antisymmetric<T,Tensor_Dim>,
+ T,Dim,i,j>(*this);
+ }
+
+ /* This is for expressions where a number is used for one slot, and
+ an index for another, yielding a Tensor1_Expr. The non-const
+ versions don't actually create a Tensor2_number_rhs_[01] object.
+ They create a Tensor1_Expr directly, which provides the
+ appropriate indexing operators. The const versions do create a
+ Tensor2_number_[01]. */
+
+ template<char i, int N, int Dim>
+ Tensor1_Expr<Tensor2_number_rhs_1<Tensor2_antisymmetric<T,Tensor_Dim>,T,N>,
+ T,Dim,i>
+ operator()(const Index<i,Dim> index1, const Number<N> &n1)
+ {
+ typedef Tensor2_number_rhs_1<Tensor2_antisymmetric<T,Tensor_Dim>,T,N>
+ TensorExpr;
+ return Tensor1_Expr<TensorExpr,T,Dim,i>(*this);
+ }
+
+ template<char i, int N, int Dim>
+ Tensor1_Expr<const Tensor2_number_1<const Tensor2_antisymmetric<T,Tensor_Dim>,
+ T,N>,T,Dim,i>
+ operator()(const Index<i,Dim> index1, const Number<N> &n1) const
+ {
+ typedef const Tensor2_number_1<const Tensor2_antisymmetric<T,Tensor_Dim>,T,N>
+ TensorExpr;
+ return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this));
+ }
+
+ template<char i, int N, int Dim>
+ Tensor1_Expr<Tensor2_number_rhs_0<Tensor2_antisymmetric<T,Tensor_Dim>,T,N>,
+ T,Dim,i>
+ operator()(const Number<N> &n1, const Index<i,Dim> index1)
+ {
+ typedef Tensor2_number_rhs_0<Tensor2_antisymmetric<T,Tensor_Dim>,T,N>
+ TensorExpr;
+ return Tensor1_Expr<TensorExpr,T,Dim,i>(*this);
+ }
+
+ template<char i, int N, int Dim>
+ Tensor1_Expr<const Tensor2_number_0<const Tensor2_antisymmetric<T,Tensor_Dim>,
+ T,N>,T,Dim,i>
+ operator()(const Number<N> &n1, const Index<i,Dim> index1) const
+ {
+ typedef const Tensor2_number_0<const Tensor2_antisymmetric<T,Tensor_Dim>,T,N>
+ TensorExpr;
+ return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this));
+ }
+
+ /* Specializations for using actual numbers instead of Number<> */
+
+ template<char i, int Dim>
+ Tensor1_Expr<const Tensor2_numeral_1<const Tensor2_antisymmetric<T,Tensor_Dim>,
+ T>,T,Dim,i>
+ operator()(const Index<i,Dim> index1, const int N) const
+ {
+ typedef const Tensor2_numeral_1<const Tensor2_antisymmetric<T,Tensor_Dim>,T>
+ TensorExpr;
+ return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this,N));
+ }
+
+ template<char i, int Dim>
+ Tensor1_Expr<const Tensor2_numeral_0<const Tensor2_antisymmetric<T,Tensor_Dim>,
+ T>,T,Dim,i>
+ operator()(const int N, const Index<i,Dim> index1) const
+ {
+ typedef const Tensor2_numeral_0<const Tensor2_antisymmetric<T,Tensor_Dim>,T>
+ TensorExpr;
+ return Tensor1_Expr<TensorExpr,T,Dim,i>(TensorExpr(*this,N));
+ }
+
+ /* These two operator()'s return the Tensor2 with internal
+ contractions, yielding a T. I have to specify one for both
+ const and non-const because otherwise the compiler will use the
+ operator() which gives a Tensor2_Expr<>. */
+
+ template<char i, int Dim>
+ T operator()(const Index<i,Dim> index1, const Index<i,Dim> index2)
+ {
+ return 0;
+ }
+
+ template<char i, int Dim>
+ T operator()(const Index<i,Dim> index1, const Index<i,Dim> index2) const
+ {
+ return 0;
+ }
+};
More information about the CIG-COMMITS
mailing list