Chaste  Release::3.4
FourthOrderTensor.hpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #ifndef _FOURTHORDERTENSOR_HPP_
37 #define _FOURTHORDERTENSOR_HPP_
38 
39 #include <cassert>
40 #include <vector>
41 
42 #include "UblasIncludes.hpp"
43 #include "Exception.hpp"
44 
50 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
52 {
53 private:
54 
55  std::vector<double> mData;
63  unsigned GetVectorIndex(unsigned M, unsigned N, unsigned P, unsigned Q)
64  {
65  assert(M<DIM1);
66  assert(N<DIM2);
67  assert(P<DIM3);
68  assert(Q<DIM4);
69  return M + DIM1*N + DIM1*DIM2*P + DIM1*DIM2*DIM3*Q;
70  }
71 
72 public:
73 
78 
87  template<unsigned CONTRACTED_DIM>
88  void SetAsContractionOnFirstDimension(const c_matrix<double,DIM1,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<CONTRACTED_DIM,DIM2,DIM3,DIM4>& rTensor);
89 
90 
99  template<unsigned CONTRACTED_DIM>
100  void SetAsContractionOnSecondDimension(const c_matrix<double,DIM2,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,CONTRACTED_DIM,DIM3,DIM4>& rTensor);
101 
110  template<unsigned CONTRACTED_DIM>
111  void SetAsContractionOnThirdDimension(const c_matrix<double,DIM3,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,DIM2,CONTRACTED_DIM,DIM4>& rTensor);
112 
121  template<unsigned CONTRACTED_DIM>
122  void SetAsContractionOnFourthDimension(const c_matrix<double,DIM4,CONTRACTED_DIM>& rMatrix, FourthOrderTensor<DIM1,DIM2,DIM3,CONTRACTED_DIM>& rTensor);
123 
132  double& operator()(unsigned M, unsigned N, unsigned P, unsigned Q);
133 
137  void Zero();
138 
142  std::vector<double>& rGetData()
143  {
144  return mData;
145  }
146 };
147 
149 // Implementation (lots of possibilities for the dimensions so no point with explicit instantiation)
151 
152 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
154 {
155  unsigned size = DIM1*DIM2*DIM3*DIM4;
156  mData.resize(size, 0.0);
157 }
158 
159 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
160 template<unsigned CONTRACTED_DIM>
162 {
163  Zero();
164 
165  std::vector<double>::iterator iter = mData.begin();
166  std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
167 
168  for (unsigned d=0; d<DIM4; d++)
169  {
170  for (unsigned c=0; c<DIM3; c++)
171  {
172  for (unsigned b=0; b<DIM2; b++)
173  {
174  for (unsigned a=0; a<DIM1; a++)
175  {
176  for (unsigned N=0; N<CONTRACTED_DIM; N++)
177  {
178  /*
179  * The following just does
180  *
181  * mData[GetVectorIndex(a,b,c,d)] += rMatrix(a,N) * rTensor(N,b,c,d);
182  *
183  * but more efficiently using iterators into the data vector, not
184  * using random access.
185  */
186  *iter += rMatrix(a,N) * *other_tensor_iter;
187  other_tensor_iter++;
188  }
189 
190  iter++;
191 
192  if (a != DIM1-1)
193  {
194  other_tensor_iter -= CONTRACTED_DIM;
195  }
196  }
197  }
198  }
199  }
200 }
201 
202 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
203 template<unsigned CONTRACTED_DIM>
205 {
206  Zero();
207 
208  std::vector<double>::iterator iter = mData.begin();
209  std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
210 
211  for (unsigned d=0; d<DIM4; d++)
212  {
213  for (unsigned c=0; c<DIM3; c++)
214  {
215  for (unsigned b=0; b<DIM2; b++)
216  {
217  for (unsigned N=0; N<CONTRACTED_DIM; N++)
218  {
219  for (unsigned a=0; a<DIM1; a++)
220  {
221  /*
222  * The following just does
223  *
224  * mData[GetVectorIndex(a,b,c,d)] += rMatrix(b,N) * rTensor(a,N,c,d);
225  *
226  * but more efficiently using iterators into the data vector, not
227  * using random access.
228  */
229  *iter += rMatrix(b,N) * *other_tensor_iter;
230  iter++;
231  other_tensor_iter++;
232  }
233 
234  if (N != CONTRACTED_DIM-1)
235  {
236  iter -= DIM1;
237  }
238  }
239  if (b != DIM2-1)
240  {
241  other_tensor_iter -= CONTRACTED_DIM*DIM1;
242  }
243  }
244  }
245  }
246 }
247 
248 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
249 template<unsigned CONTRACTED_DIM>
251 {
252  Zero();
253 
254  std::vector<double>::iterator iter = mData.begin();
255  std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
256 
257  for (unsigned d=0; d<DIM4; d++)
258  {
259  for (unsigned c=0; c<DIM3; c++)
260  {
261  for (unsigned N=0; N<CONTRACTED_DIM; N++)
262  {
263  for (unsigned b=0; b<DIM2; b++)
264  {
265  for (unsigned a=0; a<DIM1; a++)
266  {
267  /*
268  * The following just does
269  *
270  * mData[GetVectorIndex(a,b,c,d)] += rMatrix(c,N) * rTensor(a,b,N,d);
271  *
272  * but more efficiently using iterators into the data vector, not
273  * using random access.
274  */
275  *iter += rMatrix(c,N) * *other_tensor_iter;
276  iter++;
277  other_tensor_iter++;
278  }
279  }
280 
281  if (N != CONTRACTED_DIM-1)
282  {
283  iter -= DIM1*DIM2;
284  }
285  }
286 
287  if (c != DIM3-1)
288  {
289  other_tensor_iter -= CONTRACTED_DIM*DIM1*DIM2;
290  }
291  }
292  }
293 }
294 
295 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
296 template<unsigned CONTRACTED_DIM>
298 {
299  Zero();
300 
301  std::vector<double>::iterator iter = mData.begin();
302  std::vector<double>::iterator other_tensor_iter = rTensor.rGetData().begin();
303 
304  for (unsigned d=0; d<DIM4; d++)
305  {
306  for (unsigned N=0; N<CONTRACTED_DIM; N++)
307  {
308  for (unsigned c=0; c<DIM3; c++)
309  {
310  for (unsigned b=0; b<DIM2; b++)
311  {
312  for (unsigned a=0; a<DIM1; a++)
313  {
314  /*
315  * The following just does
316  *
317  * mData[GetVectorIndex(a,b,c,d)] += rMatrix(d,N) * rTensor(a,b,c,N);
318  *
319  * but more efficiently using iterators into the data vector, not
320  * using random access.
321  */
322  *iter += rMatrix(d,N) * *other_tensor_iter;
323 
324  iter++;
325  other_tensor_iter++;
326  }
327  }
328  }
329 
330  if (N != CONTRACTED_DIM-1)
331  {
332  iter-= DIM1*DIM2*DIM3;
333  }
334  }
335 
336  other_tensor_iter -= CONTRACTED_DIM*DIM1*DIM2*DIM3;
337  }
338 }
339 
340 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
341 double& FourthOrderTensor<DIM1,DIM2,DIM3,DIM4>::operator()(unsigned M, unsigned N, unsigned P, unsigned Q)
342 {
343  assert(M<DIM1);
344  assert(N<DIM2);
345  assert(P<DIM3);
346  assert(Q<DIM4);
347 
348  return mData[GetVectorIndex(M,N,P,Q)];
349 }
350 
351 template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
353 {
354  for (unsigned i=0; i<mData.size(); i++)
355  {
356  mData[i] = 0.0;
357  }
358 }
359 
360 #endif //_FOURTHORDERTENSOR_HPP_
std::vector< double > mData
double & operator()(unsigned M, unsigned N, unsigned P, unsigned Q)
std::vector< double > & rGetData()
void SetAsContractionOnFirstDimension(const c_matrix< double, DIM1, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< CONTRACTED_DIM, DIM2, DIM3, DIM4 > &rTensor)
void SetAsContractionOnFourthDimension(const c_matrix< double, DIM4, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, DIM2, DIM3, CONTRACTED_DIM > &rTensor)
void SetAsContractionOnSecondDimension(const c_matrix< double, DIM2, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, CONTRACTED_DIM, DIM3, DIM4 > &rTensor)
void SetAsContractionOnThirdDimension(const c_matrix< double, DIM3, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, DIM2, CONTRACTED_DIM, DIM4 > &rTensor)
unsigned GetVectorIndex(unsigned M, unsigned N, unsigned P, unsigned Q)