Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
FourthOrderTensor.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF 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
50template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
52{
53private:
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
72public:
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
152template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
154{
155 unsigned size = DIM1*DIM2*DIM3*DIM4;
156 mData.resize(size, 0.0);
157}
158
159template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
160template<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
202template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
203template<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
248template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
249template<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
295template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
296template<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
340template<unsigned DIM1, unsigned DIM2, unsigned DIM3, unsigned DIM4>
341double& 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
351template<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_
double & operator()(unsigned M, unsigned N, unsigned P, unsigned Q)
std::vector< double > & rGetData()
std::vector< double > mData
unsigned GetVectorIndex(unsigned M, unsigned N, unsigned P, unsigned Q)
void SetAsContractionOnSecondDimension(const c_matrix< double, DIM2, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, CONTRACTED_DIM, DIM3, DIM4 > &rTensor)
void SetAsContractionOnFourthDimension(const c_matrix< double, DIM4, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, DIM2, DIM3, CONTRACTED_DIM > &rTensor)
void SetAsContractionOnFirstDimension(const c_matrix< double, DIM1, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< CONTRACTED_DIM, DIM2, DIM3, DIM4 > &rTensor)
void SetAsContractionOnThirdDimension(const c_matrix< double, DIM3, CONTRACTED_DIM > &rMatrix, FourthOrderTensor< DIM1, DIM2, CONTRACTED_DIM, DIM4 > &rTensor)