Chaste  Release::3.4
UblasCustomFunctions.hpp
Go to the documentation of this file.
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 UBLASCUSTOMFUNCTIONS_HPP_
37 #define UBLASCUSTOMFUNCTIONS_HPP_
38 
45 #include "UblasIncludes.hpp"
46 
47 #include "Exception.hpp"
48 #include "MathsCustomFunctions.hpp"
49 
50 // COMMON DETERMINANTS - SQUARE
51 
58 template<class T>
59 inline T Determinant(const boost::numeric::ublas::c_matrix<T, 1, 1>& rM)
60 {
61  using namespace boost::numeric::ublas;
62 
63  return rM(0,0);
64 }
65 
72 template<class T>
73 T Determinant(const boost::numeric::ublas::c_matrix<T,2,2>& rM)
74 {
75  using namespace boost::numeric::ublas;
76 
77  return rM(0,0)*rM(1,1) - rM(1,0)*rM(0,1);
78 }
79 
86 template<class T>
87 T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 3>& rM)
88 {
89  using namespace boost::numeric::ublas;
90 
91  return rM(0,0) * (rM(1,1)*rM(2,2) - rM(1,2)*rM(2,1))
92  - rM(0,1) * (rM(1,0)*rM(2,2) - rM(1,2)*rM(2,0))
93  + rM(0,2) * (rM(1,0)*rM(2,1) - rM(1,1)*rM(2,0));
94 }
95 
96 // COMMON GENERALIZED DETERMINANTS - NOT SQUARE
97 
105 template<class T>
106 T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 2>& rM)
107 {
108  using namespace boost::numeric::ublas;
109  c_matrix<T,2,2> product = prod(trans(rM), rM);
110  return std::sqrt(Determinant(product));
111 }
112 
120 template<class T>
121 T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 1>& rM)
122 {
123  using namespace boost::numeric::ublas;
124  return std::sqrt(rM(0,0)*rM(0,0) + rM(1,0)*rM(1,0) + rM(2,0)*rM(2,0));
125 }
126 
134 template<class T>
135 T Determinant(const boost::numeric::ublas::c_matrix<T, 2, 1>& rM)
136 {
137  using namespace boost::numeric::ublas;
138  return std::sqrt(rM(0,0) * rM(0,0) + rM(1,0) * rM(1,0));
139 }
140 
147 template<class T>
148 T Determinant(const boost::numeric::ublas::c_matrix<T, 3, 0>& rM)
149 {
151 }
152 
159 template<class T>
160 T Determinant(const boost::numeric::ublas::c_matrix<T, 2, 0>& rM)
161 {
163 }
164 
171 template<class T>
172 T Determinant(const boost::numeric::ublas::c_matrix<T, 1, 0>& rM)
173 {
175 }
176 
177 // COMMON SUBDETERMINANTS - SQUARE
178 
188 template<class T>
189 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 1, 1>& rM, const unsigned missrow, const unsigned misscol)
190 {
191  using namespace boost::numeric::ublas;
192 
193  assert(missrow == 0);
194  assert(misscol == 0);
195  return 1.0;
196 }
197 
206 template<class T>
207 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 2, 2>& rM, const unsigned missrow, const unsigned misscol)
208 {
209  using namespace boost::numeric::ublas;
210 
211  assert(missrow < 2);
212  assert(misscol < 2);
213 
214  unsigned row = (missrow==1) ? 0 : 1;
215  unsigned col = (misscol==1) ? 0 : 1;
216  return rM(row,col);
217 }
218 
227 template<class T>
228 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 3>& rM, const unsigned missrow, const unsigned misscol)
229 {
230  using namespace boost::numeric::ublas;
231 
232  assert(missrow < 3);
233  assert(misscol < 3);
234 
235  unsigned lorow = (missrow==0) ? 1 : 0;
236  unsigned hirow = (missrow==2) ? 1 : 2;
237  unsigned locol = (misscol==0) ? 1 : 0;
238  unsigned hicol = (misscol==2) ? 1 : 2;
239  return rM(lorow,locol)*rM(hirow,hicol) - rM(lorow,hicol)*rM(hirow,locol);
240 }
241 
242 // COMMON SUBDETERMINANTS - NOT SQUARE
243 
252 template<class T>
253 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 2>& rM, const unsigned missrow, const unsigned misscol)
254 {
255  using namespace boost::numeric::ublas;
256 
257  assert(missrow < 3);
258  //assert(misscol < 2); //Don't assert this since it is used
259 
260  unsigned lorow = (missrow==0) ? 1 : 0;
261  unsigned hirow = (missrow==2) ? 1 : 2;
262  unsigned locol = (misscol==0) ? 1 : 0;
263  unsigned hicol = (misscol==2) ? 1 : 2;
264  return rM(lorow,locol)*rM(hirow,hicol) - rM(lorow,hicol)*rM(hirow,locol);
265 }
266 
275 template<class T>
276 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 1>& rM, const unsigned missrow, const unsigned misscol)
277 {
278  using namespace boost::numeric::ublas;
279 
280  assert(missrow < 3);
281  assert(misscol < 1);
282 
283  unsigned lorow = (missrow==0) ? 1 : 0;
284  unsigned hirow = (missrow==2) ? 1 : 2;
285  unsigned locol = (misscol==0) ? 1 : 0;
286  unsigned hicol = (misscol==2) ? 1 : 2;
287  return rM(lorow,locol)*rM(hirow,hicol) - rM(lorow,hicol)*rM(hirow,locol);
288 }
289 
298 template<class T>
299 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 2, 1>& rM, const unsigned missrow, const unsigned misscol)
300 {
301  using namespace boost::numeric::ublas;
302 
303  assert(missrow < 2);
304  assert(misscol < 1);
305 
306  unsigned row = (missrow==1) ? 0 : 1;
307  unsigned col = (misscol==1) ? 0 : 1;
308  return rM(row,col);
309 }
310 
311 #if defined(__xlC__)
312 /* IBM compiler doesn't support zero-sized arrays*/
313 #else //#if defined(__xlC__)
314 
322 template<class T>
323 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 3, 0>& rM, const unsigned missrow, const unsigned misscol)
324 {
326 }
327 
336 template<class T>
337 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 2, 0>& rM, const unsigned missrow, const unsigned misscol)
338 {
340 }
341 
350 template<class T>
351 T SubDeterminant(const boost::numeric::ublas::c_matrix<T, 1, 0>& rM, const unsigned missrow, const unsigned misscol)
352 {
354 }
355 #endif //#if defined(__xlC__)
356 
357 // COMMON INVERSES - SQUARE
358 
366 template<class T>
367 boost::numeric::ublas::c_matrix<T, 1, 1> Inverse(const boost::numeric::ublas::c_matrix<T, 1, 1>& rM)
368 {
369  using namespace boost::numeric::ublas;
370 
371  c_matrix<T,1,1> inverse;
372  T det = Determinant(rM);
373  assert(fabs(det) > DBL_EPSILON); // else it is a singular matrix
374  inverse(0,0) = 1.0/det;
375  return inverse;
376 }
377 
385 template<class T>
386 boost::numeric::ublas::c_matrix<T, 2, 2> Inverse(const boost::numeric::ublas::c_matrix<T, 2, 2>& rM)
387 {
388  using namespace boost::numeric::ublas;
389 
390  c_matrix<T, 2, 2> inverse;
391  T det = Determinant(rM);
392 
393  assert( fabs(det) > DBL_EPSILON ); // else it is a singular matrix
394  inverse(0,0) = rM(1,1)/det;
395  inverse(0,1) = -rM(0,1)/det;
396  inverse(1,0) = -rM(1,0)/det;
397  inverse(1,1) = rM(0,0)/det;
398  return inverse;
399 }
400 
408 template<class T>
409 boost::numeric::ublas::c_matrix<T, 3, 3> Inverse(const boost::numeric::ublas::c_matrix<T, 3, 3>& rM)
410 {
411  using namespace boost::numeric::ublas;
412 
413  c_matrix<T, 3, 3> inverse;
414  T det = Determinant(rM);
415  assert(fabs(det) > DBL_EPSILON); // else it is a singular matrix
416 
417  inverse(0,0) = (rM(1,1)*rM(2,2) - rM(1,2)*rM(2,1))/det;
418  inverse(1,0) = -(rM(1,0)*rM(2,2) - rM(1,2)*rM(2,0))/det;
419  inverse(2,0) = (rM(1,0)*rM(2,1) - rM(1,1)*rM(2,0))/det;
420  inverse(0,1) = -(rM(0,1)*rM(2,2) - rM(0,2)*rM(2,1))/det;
421  inverse(1,1) = (rM(0,0)*rM(2,2) - rM(0,2)*rM(2,0))/det;
422  inverse(2,1) = -(rM(0,0)*rM(2,1) - rM(0,1)*rM(2,0))/det;
423  inverse(0,2) = (rM(0,1)*rM(1,2) - rM(0,2)*rM(1,1))/det;
424  inverse(1,2) = -(rM(0,0)*rM(1,2) - rM(0,2)*rM(1,0))/det;
425  inverse(2,2) = (rM(0,0)*rM(1,1) - rM(0,1)*rM(1,0))/det;
426 
427  return inverse;
428 }
429 
430 // COMMON PSEUDO-INVERSES - NOT SQUARE
431 
439 template<class T>
440 boost::numeric::ublas::c_matrix<T, 2, 3> Inverse(const boost::numeric::ublas::c_matrix<T, 3, 2>& rM)
441 {
442  using namespace boost::numeric::ublas;
443 
444  c_matrix<T, 2, 3> inverse;
445 
446  //
447  // calculate (T'T)^-1, where T'T = (a b)
448  // (c d)
449 
450  T a = rM(0,0)*rM(0,0) + rM(1,0)*rM(1,0) + rM(2,0)*rM(2,0);
451  T b = rM(0,0)*rM(0,1) + rM(1,0)*rM(1,1) + rM(2,0)*rM(2,1);
452  T c = b;
453  T d = rM(0,1)*rM(0,1) + rM(1,1)*rM(1,1) + rM(2,1)*rM(2,1);
454 
455  T det = a*d - b*c;
456 
457  T a_inv = d/det;
458  T b_inv = -b/det;
459  T c_inv = -c/det;
460  T d_inv = a/det;
461 
462  inverse(0,0) = a_inv*rM(0,0) + b_inv*rM(0,1);
463  inverse(1,0) = c_inv*rM(0,0) + d_inv*rM(0,1);
464  inverse(0,1) = a_inv*rM(1,0) + b_inv*rM(1,1);
465  inverse(1,1) = c_inv*rM(1,0) + d_inv*rM(1,1);
466  inverse(0,2) = a_inv*rM(2,0) + b_inv*rM(2,1);
467  inverse(1,2) = c_inv*rM(2,0) + d_inv*rM(2,1);
468 
469  return inverse;
470 }
471 
479 template<class T>
480 boost::numeric::ublas::c_matrix<T, 1, 2> Inverse(const boost::numeric::ublas::c_matrix<T, 2, 1>& rM)
481 {
482  using namespace boost::numeric::ublas;
483 
484  c_matrix<T, 1, 2> inverse;
485  T det = Determinant(rM);
486 
487  inverse(0,0) = rM(0,0)/det/det;
488  inverse(0,1) = rM(1,0)/det/det;
489 
490  return inverse;
491 }
492 
500 template<class T>
501 boost::numeric::ublas::c_matrix<T, 1, 3> Inverse(const boost::numeric::ublas::c_matrix<T, 3, 1>& rM)
502 {
503  using namespace boost::numeric::ublas;
504 
505  c_matrix<T, 1, 3> inverse;
506  T det = Determinant(rM);
507 
508  inverse(0,0) = rM(0,0)/det/det;
509  inverse(0,1) = rM(1,0)/det/det;
510  inverse(0,2) = rM(2,0)/det/det;
511 
512  return inverse;
513 }
514 
515 // COMMON MATRIX TRACES
516 
523 template<class T>
524 inline T Trace(const c_matrix<T, 1, 1>& rM)
525 {
526  return rM(0,0);
527 }
528 
535 template<class T>
536 inline T Trace(const c_matrix<T, 2, 2>& rM)
537 {
538  return rM(0,0) + rM(1,1);
539 }
540 
547 template<class T>
548 inline T Trace(const c_matrix<T, 3, 3>& rM)
549 {
550  return rM(0,0) + rM(1,1) + rM(2,2);
551 }
552 
559 template<class T>
560 inline T Trace(const c_matrix<T, 4, 4>& rM)
561 {
562  return rM(0,0) + rM(1,1) + rM(2,2) + rM(3,3);
563 }
564 
565 // OTHER MATRIX FUNCTIONS (INVARIANTS, EIGENVECTORS)
566 
575 template<class T>
576 T SecondInvariant(const c_matrix<T, 3, 3>& rM)
577 {
578  return rM(0,0)*rM(1,1) + rM(1,1)*rM(2,2) + rM(2,2)*rM(0,0)
579  - rM(1,0)*rM(1,0) - rM(2,1)*rM(2,1) - rM(2,0)*rM(2,0);
580 }
581 
589 template<class T>
590 T SecondInvariant(const c_matrix<T, 2, 2>& rM)
591 {
592  return Determinant(rM);
593 }
594 
604 c_vector<double,3> CalculateEigenvectorForSmallestNonzeroEigenvalue(c_matrix<double, 3, 3>& rA);
605 
612 double CalculateMaxEigenpair(c_matrix<double, 3, 3>& rA, c_vector<double, 3>& rEigenvector);
613 
614  //COMMON VECTOR FUNCTIONS
615 
616 
625 template<class T>
626 c_vector<T, 1> VectorProduct(const c_vector<T, 1>& rA, const c_vector<T, 1>& rB)
627 {
629 }
638 template<class T>
639 c_vector<T, 2> VectorProduct(const c_vector<T, 2>& rA, const c_vector<T, 2>& rB)
640 {
642 }
643 
651 template<class T>
652 c_vector<T, 3> VectorProduct(const c_vector<T, 3>& rA, const c_vector<T, 3>& rB)
653 {
654 
655  c_vector<T, 3> result;
656 
657  result(0) = rA(1)*rB(2) - rA(2)*rB(1);
658  result(1) = rA(2)*rB(0) - rA(0)*rB(2);
659  result(2) = rA(0)*rB(1) - rA(1)*rB(0);
660 
661  return result;
662 }
663 
670 c_vector<double, 1> Create_c_vector(double x);
671 
679 c_vector<double, 2> Create_c_vector(double x, double y);
680 
689 c_vector<double, 3> Create_c_vector(double x, double y, double z);
690 
691 #endif /*UBLASCUSTOMFUNCTIONS_HPP_*/
T SecondInvariant(const c_matrix< T, 3, 3 > &rM)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
c_vector< double, 3 > CalculateEigenvectorForSmallestNonzeroEigenvalue(c_matrix< double, 3, 3 > &rA)
T Trace(const c_matrix< T, 1, 1 > &rM)
c_vector< T, 1 > VectorProduct(const c_vector< T, 1 > &rA, const c_vector< T, 1 > &rB)
T Determinant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
#define NEVER_REACHED
Definition: Exception.hpp:206
c_vector< double, 1 > Create_c_vector(double x)
T SubDeterminant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM, const unsigned missrow, const unsigned misscol)
double CalculateMaxEigenpair(c_matrix< double, 3, 3 > &rA, c_vector< double, 3 > &rEigenvector)