tensor4.cc
Go to the documentation of this file.
#include "rheolef/tensor4.h"
namespace rheolef {
// identity tensor:
// tensor y = ddot(Id,x) ==> x==y
template<class T>
tensor4_basic<T>
tensor4_basic<T>::eye (size_type d)
{
tensor4_basic<T> a;
for (size_t i = 0; i < d; i++)
for (size_t j = 0; j < d; j++) a(i,j,i,j) = 1;
return a;
}
// algebra
template<class T>
tensor_basic<T>
ddot (const tensor4_basic<T>& a, const tensor_basic<T>& x)
{
typedef typename tensor4_basic<T>::size_type size_type;
tensor_basic<T> y;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
y(i,j) += a(i,j,k,l)*x(k,l);
return y;
}
template<class T>
tensor_basic<T>
ddot (const tensor_basic<T>& y, const tensor4_basic<T>& a)
{
typedef typename tensor4_basic<T>::size_type size_type;
tensor_basic<T> x;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
x(k,l) += y(i,j)*a(i,j,k,l);
return x;
}
template<class T>
tensor4_basic<T>::tensor4_basic (
const std::initializer_list<std::initializer_list<
std::initializer_list<std::initializer_list<T> > > >& il)
: _x (tensor_basic<T>(T()))
{
typedef std::initializer_list<T> L1;
typedef std::initializer_list<L1> L2;
typedef std::initializer_list<L2> L3;
typedef std::initializer_list<L3> L4;
typedef typename std::initializer_list<L4>::const_iterator const_iterator;
typedef typename std::initializer_list<L3>::const_iterator const_iterator_row3;
typedef typename std::initializer_list<L2>::const_iterator const_iterator_row2;
typedef typename std::initializer_list<L1>::const_iterator const_iterator_row1;
typedef typename std::initializer_list<T>::const_iterator const_iterator_elt;
this->operator= (T());
check_macro (il.size() <= 3, "unexpected initializer list size=" << il.size() << " > 3");
size_type i = 0;
for (const_iterator_row3 iter = il.begin(); iter != il.end(); ++iter, ++i) {
const L3& row3 = *iter;
check_macro (row3.size() <= 3, "unexpected initializer list size=" << row3.size() << " > 3");
size_type j = 0;
for (const_iterator_row2 jter = row3.begin(); jter != row3.end(); ++jter, ++j) {
const L2& row2 = *jter;
check_macro (row2.size() <= 3, "unexpected initializer list size=" << row2.size() << " > 3");
size_type k = 0;
for (const_iterator_row1 kter = row2.begin(); kter != row2.end(); ++kter, ++k) {
const L1& row1 = *kter;
check_macro (row1.size() <= 3, "unexpected initializer list size=" << row1.size() << " > 3");
size_type l = 0;
for (const_iterator_elt lter = row1.begin(); lter != row1.end(); ++lter, ++l) {
const T& elt = *lter;
operator() (i,j,k,l) = elt;
}
}
}
}
}
template<class T>
tensor4_basic<T>&
tensor4_basic<T>::operator= (const T& val)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
operator() (i,j,k,l) = val;
return *this;
}
template<class T>
tensor4_basic<T>&
tensor4_basic<T>::operator= (const tensor4_basic<T>& a)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
operator() (i,j,k,l) = a(i,j,k,l);
return *this;
}
template<class T>
tensor4_basic<T>
tensor4_basic<T>::operator+ (const tensor4_basic<T>& b) const
{
typedef typename tensor4_basic<T>::size_type size_type;
tensor4_basic<T> c;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
c(i,j,k,l) = operator() (i,j,k,l) + b(i,j,k,l);
return c;
}
template<class T>
tensor4_basic<T>
tensor4_basic<T>::operator- (const tensor4_basic<T>& b) const
{
typedef typename tensor4_basic<T>::size_type size_type;
tensor4_basic<T> c;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
c(i,j,k,l) = operator() (i,j,k,l) - b(i,j,k,l);
return c;
}
template<class T>
tensor4_basic<T>&
tensor4_basic<T>::operator+= (const tensor4_basic<T>& b)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
operator() (i,j,k,l) += b(i,j,k,l);
return *this;
}
template<class T>
tensor4_basic<T>&
tensor4_basic<T>::operator-= (const tensor4_basic<T>& b)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
operator() (i,j,k,l) -= b(i,j,k,l);
return *this;
}
template<class T>
tensor4_basic<T>&
tensor4_basic<T>::operator*= (const T& fact)
{
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
operator() (i,j,k,l) *= fact;
return *this;
}
template <class T>
T
norm2 (const tensor4_basic<T>& a)
{
typedef typename tensor4_basic<T>::size_type size_type;
T sum = 0;
for (size_type i = 0; i < 3; i++)
for (size_type j = 0; j < 3; j++)
for (size_type k = 0; k < 3; k++)
for (size_type l = 0; l < 3; l++)
sum += sqr(a(i,j,k,l));
return sum;
}
// output
template<class T>
std::ostream&
tensor4_basic<T>::put (std::ostream& out, size_type d) const
{
using namespace std;
switch (d) {
case 0 : return out;
case 1 : return out << _x(0,0)(0,0);
case 2 : return out << "[[" << _x(0,0)(0,0) << ", " << _x(0,0)(0,1) << ";" << endl
<< " " << _x(0,0)(1,0) << ", " << _x(0,0)(1,1) << "]," << endl
<< " [" << _x(0,1)(0,0) << ", " << _x(0,1)(0,1) << ";" << endl
<< " " << _x(0,1)(1,0) << ", " << _x(0,1)(1,1) << "];" << endl
<< " [" << _x(1,0)(0,0) << ", " << _x(1,0)(0,1) << ";" << endl
<< " " << _x(1,0)(1,0) << ", " << _x(1,0)(1,1) << "]," << endl
<< " [" << _x(1,1)(0,0) << ", " << _x(1,1)(0,1) << ";" << endl
<< " " << _x(1,1)(1,0) << ", " << _x(1,1)(1,1) << "]]" << endl;
default: fatal_macro("put(d="<<d<<"): not yet, sorry"); return out;
}
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class tensor4_basic<T>; \
template tensor_basic<T> ddot (const tensor4_basic<T>&, const tensor_basic<T>&); \
template tensor_basic<T> ddot (const tensor_basic<T>&, const tensor4_basic<T>&); \
template T norm2 (const tensor4_basic<T>&);
_RHEOLEF_instanciation(Float)
}// namespace rheolef
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
rheolef::std
Definition: vec_expr_v2.h:402
rheolef::tensor4_basic
Definition: tensor4.h:80
rheolef::tensor4_basic::put
std::ostream & put(std::ostream &out, size_type d=3) const
Definition: tensor4.cc:198
rheolef::tensor4_basic::operator+
tensor4_basic< T > operator+(const tensor4_basic< T > &b) const
Definition: tensor4.cc:125
rheolef::tensor4_basic::operator-
tensor4_basic< T > operator-(const tensor4_basic< T > &b) const
Definition: tensor4.cc:138
rheolef::tensor4_basic::size_type
size_t size_type
Definition: tensor4.h:83
rheolef::tensor4_basic::operator-=
tensor4_basic< T > & operator-=(const tensor4_basic< T > &)
Definition: tensor4.cc:162
rheolef::tensor4_basic::tensor4_basic
tensor4_basic()
Definition: tensor4.h:160
rheolef::tensor4_basic::operator+=
tensor4_basic< T > & operator+=(const tensor4_basic< T > &)
Definition: tensor4.cc:151
rheolef::tensor4_basic::operator*=
tensor4_basic< T > & operator*=(const T &k)
Definition: tensor4.cc:173
rheolef::tensor4_basic::operator=
tensor4_basic< T > & operator=(const tensor4_basic< T > &a)
Definition: tensor4.cc:114
rheolef::tensor4_basic::eye
static tensor4_basic< T > eye(size_type d=3)
Definition: tensor4.cc:29
rheolef::tensor4_basic::operator()
T & operator()(size_type i, size_type j, size_type k, size_type l)
Definition: tensor4.h:194
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
rheolef::io::out
@ out
Definition: rheostream.h:167
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
rheolef::ddot
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
a
Definition: diffusion_isotropic.h:25
rheolef::const_iterator
Definition: field_expr_recursive.h:552