Chombo + EB  3.0
RealTensor.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 #ifndef _REALTENSOR_H_
12 #define _REALTENSOR_H_
13 
14 #include <algorithm>
15 
16 #include "SPACE.H"
17 #include "REAL.H"
18 #include "MayDay.H"
19 
20 #include "NamespaceHeader.H"
21 
22 //! \class RealTensor
23 //! This class stores tensors that represent outer products of RealVects.
25 {
26  public:
27 
29  {
30  std::fill(m_T, m_T + SpaceDim*SpaceDim, 0.0);
31  }
32 
33 #if CH_SPACEDIM == 1
34  explicit RealTensor(Real a_t11)
35  {
36  m_T[0] = a_t11;
37  }
38 #elif CH_SPACEDIM == 2
39  RealTensor(Real a_t11, Real a_t12,
40  Real a_t21, Real a_t22)
41  {
42  m_T[0] = a_t11; m_T[1] = a_t12;
43  m_T[2] = a_t21; m_T[3] = a_t22;
44  }
45 #elif CH_SPACEDIM == 3
46  RealTensor(Real a_t11, Real a_t12, Real a_t13,
47  Real a_t21, Real a_t22, Real a_t23,
48  Real a_t31, Real a_t32, Real a_t33)
49  {
50  m_T[0] = a_t11; m_T[1] = a_t12; m_T[2] = a_t13;
51  m_T[3] = a_t21; m_T[4] = a_t22; m_T[5] = a_t23;
52  m_T[6] = a_t31; m_T[7] = a_t32; m_T[8] = a_t33;
53  }
54 #elif CH_SPACEDIM == 4
55  RealTensor(Real a_t11, Real a_t12, Real a_t13, Real a_t14,
56  Real a_t21, Real a_t22, Real a_t23, Real a_t24,
57  Real a_t31, Real a_t32, Real a_t33, Real a_t34,
58  Real a_t41, Real a_t42, Real a_t43, Real a_t44)
59  {
60  m_T[ 0] = a_t11; m_T[ 1] = a_t12; m_T[ 2] = a_t13; m_T[ 3] = a_t14;
61  m_T[ 4] = a_t21; m_T[ 5] = a_t22; m_T[ 6] = a_t23; m_T[ 7] = a_t24;
62  m_T[ 8] = a_t31; m_T[ 9] = a_t32; m_T[10] = a_t33; m_T[11] = a_t34;
63  m_T[12] = a_t31; m_T[13] = a_t32; m_T[14] = a_t33; m_T[15] = a_t44;
64  }
65 #elif CH_SPACEDIM == 5
66  RealTensor(Real a_t11, Real a_t12, Real a_t13, Real a_t14, Real a_t15,
67  Real a_t21, Real a_t22, Real a_t23, Real a_t24, Real a_t25,
68  Real a_t31, Real a_t32, Real a_t33, Real a_t34, Real a_t35,
69  Real a_t41, Real a_t42, Real a_t43, Real a_t44, Real a_t45,
70  Real a_t51, Real a_t52, Real a_t53, Real a_t54, Real a_t55)
71  {
72  m_T[ 0] = a_t11; m_T[ 1] = a_t12; m_T[ 2] = a_t13; m_T[ 3] = a_t14; m_T[ 4] = a_t15;
73  m_T[ 5] = a_t21; m_T[ 6] = a_t22; m_T[ 7] = a_t23; m_T[ 8] = a_t24; m_T[ 9] = a_t25;
74  m_T[10] = a_t31; m_T[11] = a_t32; m_T[12] = a_t33; m_T[13] = a_t34; m_T[14] = a_t35;
75  m_T[15] = a_t41; m_T[16] = a_t42; m_T[17] = a_t43; m_T[18] = a_t44; m_T[19] = a_t45;
76  m_T[20] = a_t51; m_T[21] = a_t52; m_T[22] = a_t53; m_T[23] = a_t54; m_T[24] = a_t55;
77  }
78 #elif CH_SPACEDIM == 6
79  RealTensor(Real a_t11, Real a_t12, Real a_t13, Real a_t14, Real a_t15, Real a_t16,
80  Real a_t21, Real a_t22, Real a_t23, Real a_t24, Real a_t25, Real a_t26,
81  Real a_t31, Real a_t32, Real a_t33, Real a_t34, Real a_t35, Real a_t36,
82  Real a_t41, Real a_t42, Real a_t43, Real a_t44, Real a_t45, Real a_t46,
83  Real a_t51, Real a_t52, Real a_t53, Real a_t54, Real a_t55, Real a_t56,
84  Real a_t61, Real a_t62, Real a_t63, Real a_t64, Real a_t65, Real a_t66)
85  {
86  m_T[ 0] = a_t11; m_T[ 1] = a_t12; m_T[ 2] = a_t13; m_T[ 3] = a_t14; m_T[ 4] = a_t15; m_T[ 5] = a_t16;
87  m_T[ 6] = a_t21; m_T[ 7] = a_t22; m_T[ 8] = a_t23; m_T[ 9] = a_t24; m_T[10] = a_t25; m_T[11] = a_t26;
88  m_T[12] = a_t31; m_T[13] = a_t32; m_T[14] = a_t33; m_T[15] = a_t34; m_T[16] = a_t35; m_T[17] = a_t36;
89  m_T[18] = a_t41; m_T[19] = a_t42; m_T[20] = a_t43; m_T[21] = a_t44; m_T[22] = a_t45; m_T[23] = a_t46;
90  m_T[24] = a_t51; m_T[25] = a_t52; m_T[26] = a_t53; m_T[27] = a_t54; m_T[28] = a_t55; m_T[29] = a_t56;
91  m_T[30] = a_t61; m_T[31] = a_t62; m_T[32] = a_t63; m_T[33] = a_t64; m_T[34] = a_t65; m_T[35] = a_t66;
92  }
93 #else
94 #error "RealTensor is only implemented for SpaceDim <= 6."
95 #endif
96 
97  // Copy constructor.
98  RealTensor(const RealTensor& a_rhs)
99  {
100  std::copy(a_rhs.m_T, a_rhs.m_T + SpaceDim*SpaceDim, m_T);
101  }
102 
103  //! Destructor.
105  {
106  }
107 
108  //! Assignment operator.
110  {
111  if (this != &a_rhs)
112  {
113  std::copy(a_rhs.m_T, a_rhs.m_T + SpaceDim*SpaceDim, m_T);
114  }
115 
116  return *this;
117  }
118 
119  //! Returns the (i,j) component of this tensor.
120  Real operator()(int a_i, int a_j) const
121  {
122  return m_T[SpaceDim*a_j + a_i];
123  }
124 
125  Real& operator()(int a_i, int a_j)
126  {
127  return m_T[SpaceDim*a_j + a_i];
128  }
129 
130  //! Computes the determinant of this tensor using cofactors.
131  Real det() const
132  {
133  const Real* tensor = m_T;
134 
135  return recurDet(tensor,SpaceDim);
136  }
137 
138  //! Computes the transpose of this tensor.
140  {
141  RealTensor t = *this;
142 #if CH_SPACEDIM >= 2
143  std::swap(t(0, 1), t(1, 0));
144 #endif
145 #if CH_SPACEDIM >= 3
146  std::swap(t(0, 2), t(2, 0));
147  std::swap(t(1, 2), t(2, 1));
148 #endif
149 #if CH_SPACEDIM >= 4
150  std::swap(t(0, 3), t(3, 0));
151  std::swap(t(1, 3), t(3, 1));
152  std::swap(t(2, 3), t(3, 2));
153 #endif
154 #if CH_SPACEDIM >= 5
155  std::swap(t(0, 4), t(4, 0));
156  std::swap(t(1, 4), t(4, 1));
157  std::swap(t(2, 4), t(4, 2));
158  std::swap(t(3, 4), t(4, 3));
159 #endif
160 #if CH_SPACEDIM == 6
161  std::swap(t(0, 5), t(5, 0));
162  std::swap(t(1, 5), t(5, 1));
163  std::swap(t(2, 5), t(5, 2));
164  std::swap(t(3, 5), t(5, 3));
165  std::swap(t(4, 5), t(5, 4));
166 #endif
167  return t;
168  }
169 
170  private:
171 
172  // Compute the determinant of this tensor using cofactors.
173  Real recurDet(const Real* a_tensor,
174  const int& a_size) const
175  {
176  Real retval = 0.0;
177 
178  if (a_size == 1)
179  {
180  retval += *a_tensor;
181  }
182  else
183  {
184  int subSize = a_size - 1;
185  Real* subTensor = new Real[subSize*subSize];
186 
187  Real sign = 1.0;
188 
189  int row1 = 0;
190  for (int col1 = 0; col1 < a_size; col1++)
191  {
192  int subRow = 0;
193  for (int row2 = 0; row2 < a_size; row2++)
194  {
195  if (row2 != row1)
196  {
197  int subCol = 0;
198  for (int col2 = 0; col2 < a_size; col2++)
199  {
200  if (col2 != col1)
201  {
202  subTensor[subSize*subRow + subCol] = a_tensor[a_size*row2 + col2];
203 
204  subCol++;
205  }
206  }
207 
208  subRow++;
209  }
210  }
211 
212  retval += sign * a_tensor[a_size*row1 + col1] * recurDet(subTensor,subSize);
213 
214  sign *= -1.0;
215  }
216 
217  delete [] subTensor;
218  }
219 
220  return retval;
221  }
222 
224 };
225 
226 #include "NamespaceFooter.H"
227 
228 #endif
Definition: RealTensor.H:24
#define CH_SPACEDIM
Definition: SPACE.H:52
RealTensor()
Definition: RealTensor.H:28
Real recurDet(const Real *a_tensor, const int &a_size) const
Definition: RealTensor.H:173
const int SpaceDim
Definition: SPACE.H:39
Real det() const
Computes the determinant of this tensor using cofactors.
Definition: RealTensor.H:131
int sign(const Side::LoHiSide &a_side)
double Real
Definition: REAL.H:33
RealTensor transpose() const
Computes the transpose of this tensor.
Definition: RealTensor.H:139
Real operator()(int a_i, int a_j) const
Returns the (i,j) component of this tensor.
Definition: RealTensor.H:120
~RealTensor()
Destructor.
Definition: RealTensor.H:104
RealTensor & operator=(const RealTensor &a_rhs)
Assignment operator.
Definition: RealTensor.H:109
RealTensor(const RealTensor &a_rhs)
Definition: RealTensor.H:98
Real & operator()(int a_i, int a_j)
Definition: RealTensor.H:125
Real m_T[CH_SPACEDIM *CH_SPACEDIM]
Definition: RealTensor.H:223