1#ifndef SET_MATRIX4_MAJORMINOR_H
2#define SET_MATRIX4_MAJORMINOR_H
31 if (
uid==0 )
return data[0];
32 else if (
uid==8 ||
uid==4 ||
uid==2 ||
uid==1 )
return data[1];
33 else if (
uid==12 ||
uid==3 )
return data[2];
34 else if (
uid==10 ||
uid==6 ||
uid==9 ||
uid==5 )
return data[3];
35 else if (
uid==14 ||
uid==13 ||
uid==11 ||
uid==7 )
return data[4];
36 else if (
uid==15 )
return data[5];
44 if (
uid==0 )
return data[0];
45 else if (
uid==8 ||
uid==4 ||
uid==2 ||
uid==1 )
return data[1];
46 else if (
uid==12 ||
uid==3 )
return data[2];
47 else if (
uid==10 ||
uid==6 ||
uid==9 ||
uid==5 )
return data[3];
48 else if (
uid==14 ||
uid==13 ||
uid==11 ||
uid==7 )
return data[4];
49 else if (
uid==15 )
return data[5];
56 for (
int k = 0;
k < 2;
k++)
58 for (
int i = -1;
i < 3;
i++)
60 for (
int l = 0;
l < 2;
l++)
62 if (
i==-1)
os <<
"┌ ┐ ";
63 else if (
i == 2)
os <<
"└ ┘ ";
67 for (
int j = 0;
j < 2;
j++)
70 os << std::scientific << std::setw(11) << val ;
83 return std::sqrt(data[0]*data[0] + data[1]*data[1] + data[2]*data[2] + data[3]*data[3] + data[4]*data[4] + data[5]*data[5]);
87 if (std::isnan(data[ 0]))
return true;
88 if (std::isnan(data[ 1]))
return true;
89 if (std::isnan(data[ 2]))
return true;
90 if (std::isnan(data[ 3]))
return true;
91 if (std::isnan(data[ 4]))
return true;
92 if (std::isnan(data[ 5]))
return true;
126 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
131 for(
int i = 0;
i < 3;
i++)
132 for(
int j = 0;
j < 3;
j++)
133 for(
int k = 0;
k < 3;
k++)
134 for(
int l = 0;
l < 3;
l++)
141 for(
int p = 0;
p < 2;
p++)
142 for(
int q = 0; q < 2; q++)
143 for(
int s = 0;
s < 2;
s++)
144 for(
int t = 0;
t < 2;
t++)
147 for(
int i = 0;
i < 3;
i++)
148 for(
int j = 0;
j < 3;
j++)
149 for(
int k = 0;
k < 3;
k++)
150 for(
int l = 0;
l < 3;
l++)
159 R = Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
160 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
161 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
167 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
168 return Cubic(C11,C12,C44,R);
177 for(
int i = 0;
i < 3;
i++)
178 for(
int j = 0;
j < 3;
j++)
179 for(
int k = 0;
k < 3;
k++)
180 for(
int l = 0;
l < 3;
l++)
187 else if(
i==
j &&
k==
l)
192 else if((
i==
k &&
j==
l) && (
i!=
j))
203 Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
204 Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
206 for(
int p = 0;
p < 2;
p++)
207 for(
int q = 0; q < 2; q++)
208 for(
int s = 0;
s < 2;
s++)
209 for(
int t = 0;
t < 2;
t++)
212 for(
int i = 0;
i < 3;
i++)
213 for(
int j = 0;
j < 3;
j++)
214 for(
int k = 0;
k < 3;
k++)
215 for(
int l = 0;
l < 3;
l++)
222 Eigen::Quaterniond q =
223 Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
224 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
225 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
226 Eigen::Matrix3d R = q.toRotationMatrix();
231 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
232 return Transverse(C11,C12,
C13,
C33,C44,R);
248 for (
int i = 0 ;
i < 6 ;
i++)
if (a.data[
i] !=
b.data[
i])
return false;
255 for (
int i = 0 ;
i < 6 ;
i++)
ret.data[
i] = a.data[
i] +
b.data[
i];
262 for (
int i = 0 ;
i < 6 ;
i++)
ret.data[
i] = a.data[
i] -
b.data[
i];
269 for (
int i = 0 ;
i < 6 ;
i++)
ret.data[
i] = a.data[
i] *
b;
281 for (
int i = 0 ;
i < 6 ;
i++)
ret.data[
i] = a.data[
i] /
b;
287 Eigen::Matrix<Set::Scalar,2,2>
ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
288 ret(0,0) = a.data[0]*
b(0,0) + a.data[1]*(
b(0,1) +
b(1,0)) + a.data[2]*
b(1,1);
289 ret(0,1) = a.data[1]*
b(0,0) + a.data[3]*(
b(0,1) +
b(1,0)) + a.data[4]*
b(1,1);
290 ret(1,1) = a.data[2]*
b(0,0) + a.data[4]*(
b(0,1) +
b(1,0)) + a.data[5]*
b(1,1);
301 ret(0) = a.data[0]*
b(0,0,0) + a.data[1]*(
b(0,1,0) +
b(1,0,0) +
b(0,0,1)) + a.data[2]*
b(1,1,0) + a.data[3]*(
b(0,1,1) +
b(1,0,1)) + a.data[4]*
b(1,1,1);
302 ret(1) = a.data[1]*
b(0,0,0) + a.data[2]*
b(0,0,1) + a.data[3]*(
b(0,1,0) +
b(1,0,0)) + a.data[4]*(
b(1,1,0) +
b(0,1,1) +
b(1,0,1)) + a.data[5]*
b(1,1,1);
330 int uid =
i + 3*
j + 9*
k + 27*
l;
332 if (
uid==0 )
return data[0];
334 else if (
uid==27 ||
uid==9 ||
uid==3 ||
uid==1 )
return data[1];
336 else if (
uid==54 ||
uid==18 ||
uid==6 ||
uid==2 )
return data[2];
338 else if (
uid==36 ||
uid==4 )
return data[3];
340 else if (
uid==63 ||
uid==45 ||
uid==7 ||
uid==5 )
return data[4];
342 else if (
uid==72 ||
uid==8 )
return data[5];
344 else if (
uid==30 ||
uid==12 ||
uid==28 ||
uid==10 )
return data[6];
348 else if (
uid==39 ||
uid==37 ||
uid==31 ||
uid==13 )
return data[8];
352 else if (
uid==75 ||
uid==73 ||
uid==35 ||
uid==17 )
return data[10];
354 else if (
uid==60 ||
uid==24 ||
uid==56 ||
uid==20 )
return data[11];
356 else if (
uid==42 ||
uid==58 ||
uid==22 ||
uid==38 )
return data[12];
360 else if (
uid==78 ||
uid==74 ||
uid==62 ||
uid==26 )
return data[14];
362 else if (
uid==40 )
return data[15];
364 else if (
uid==67 ||
uid==49 ||
uid==43 ||
uid==41 )
return data[16];
366 else if (
uid==76 ||
uid==44 )
return data[17];
368 else if (
uid==70 ||
uid==52 ||
uid==68 ||
uid==50 )
return data[18];
370 else if (
uid==79 ||
uid==77 ||
uid==71 ||
uid==53 )
return data[19];
372 else if (
uid==80 )
return data[20];
379 int uid =
i + 3*
j + 9*
k + 27*
l;
381 if (
uid==0 )
return data[0];
383 else if (
uid==27 ||
uid==9 ||
uid==3 ||
uid==1 )
return data[1];
385 else if (
uid==54 ||
uid==18 ||
uid==6 ||
uid==2 )
return data[2];
387 else if (
uid==36 ||
uid==4 )
return data[3];
389 else if (
uid==63 ||
uid==45 ||
uid==7 ||
uid==5 )
return data[4];
391 else if (
uid==72 ||
uid==8 )
return data[5];
393 else if (
uid==30 ||
uid==12 ||
uid==28 ||
uid==10 )
return data[6];
397 else if (
uid==39 ||
uid==37 ||
uid==31 ||
uid==13 )
return data[8];
401 else if (
uid==75 ||
uid==73 ||
uid==35 ||
uid==17 )
return data[10];
403 else if (
uid==60 ||
uid==24 ||
uid==56 ||
uid==20 )
return data[11];
405 else if (
uid==42 ||
uid==58 ||
uid==22 ||
uid==38 )
return data[12];
409 else if (
uid==78 ||
uid==74 ||
uid==62 ||
uid==26 )
return data[14];
411 else if (
uid==40 )
return data[15];
413 else if (
uid==67 ||
uid==49 ||
uid==43 ||
uid==41 )
return data[16];
415 else if (
uid==76 ||
uid==44 )
return data[17];
417 else if (
uid==70 ||
uid==52 ||
uid==68 ||
uid==50 )
return data[18];
419 else if (
uid==79 ||
uid==77 ||
uid==71 ||
uid==53 )
return data[19];
421 else if (
uid==80 )
return data[20];
427 for (
int i = 0;
i < 21;
i++)
428 os <<
"i = " <<
i <<
" " << data[
i] << std::endl;
431 for (
int k = 0;
k < 3;
k++)
433 for (
int i = -1;
i < 4;
i++)
435 for (
int l = 0;
l < 3;
l++)
437 if (
i==-1)
os <<
"┌ ┐ ";
438 else if (
i == 3)
os <<
"└ ┘ ";
442 for (
int j = 0;
j < 3;
j++)
445 os << std::scientific << std::setw(11) << val ;
459 for (
int i = 0;
i < 21;
i++)
retsq += data[
i];
460 return std::sqrt(
retsq);
464 if (std::isnan(data[ 0]))
return true;
465 if (std::isnan(data[ 1]))
return true;
466 if (std::isnan(data[ 2]))
return true;
467 if (std::isnan(data[ 3]))
return true;
468 if (std::isnan(data[ 4]))
return true;
469 if (std::isnan(data[ 5]))
return true;
470 if (std::isnan(data[ 6]))
return true;
471 if (std::isnan(data[ 7]))
return true;
472 if (std::isnan(data[ 8]))
return true;
473 if (std::isnan(data[ 9]))
return true;
474 if (std::isnan(data[10]))
return true;
475 if (std::isnan(data[11]))
return true;
476 if (std::isnan(data[12]))
return true;
477 if (std::isnan(data[13]))
return true;
478 if (std::isnan(data[14]))
return true;
479 if (std::isnan(data[15]))
return true;
480 if (std::isnan(data[16]))
return true;
481 if (std::isnan(data[17]))
return true;
482 if (std::isnan(data[18]))
return true;
483 if (std::isnan(data[19]))
return true;
484 if (std::isnan(data[20]))
return true;
519 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
524 for(
int i = 0;
i < 3;
i++)
525 for(
int j = 0;
j < 3;
j++)
526 for(
int k = 0;
k < 3;
k++)
527 for(
int l = 0;
l < 3;
l++)
534 for(
int p = 0;
p < 3;
p++)
535 for(
int q = 0; q < 3; q++)
536 for(
int s = 0;
s < 3;
s++)
537 for(
int t = 0;
t < 3;
t++)
540 for(
int i = 0;
i < 3;
i++)
541 for(
int j = 0;
j < 3;
j++)
542 for(
int k = 0;
k < 3;
k++)
543 for(
int l = 0;
l < 3;
l++)
552 R = Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
553 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
554 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
560 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
561 return Cubic(C11,C12,C44,R);
569 for(
int i = 0;
i < 3;
i++)
570 for(
int j = 0;
j < 3;
j++)
571 for(
int k = 0;
k < 3;
k++)
572 for(
int l = 0;
l < 3;
l++)
579 else if(
i==
j &&
k==
l)
584 else if((
i==
k &&
j==
l) && (
i!=
j))
595 Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
596 Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
598 for(
int p = 0;
p < 3;
p++)
599 for(
int q = 0; q < 3; q++)
600 for(
int s = 0;
s < 3;
s++)
601 for(
int t = 0;
t < 3;
t++)
604 for(
int i = 0;
i < 3;
i++)
605 for(
int j = 0;
j < 3;
j++)
606 for(
int k = 0;
k < 3;
k++)
607 for(
int l = 0;
l < 3;
l++)
614 Eigen::Quaterniond q =
615 Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
616 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
617 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
618 Eigen::Matrix3d R = q.toRotationMatrix();
623 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
624 return Transverse(C11,C12,
C13,
C33,C44,R);
639 for (
int i = 0 ;
i < 21 ;
i++)
if (a.data[
i] !=
b.data[
i])
return false;
646 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] +
b.data[
i];
653 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] -
b.data[
i];
660 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] *
b;
672 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] /
b;
678 Eigen::Matrix<amrex::Real,3,3>
ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
679 ret(0,0) += a.data[ 0]*
b(0,0) + a.data[ 1]*(
b(0,1) +
b(1,0)) + a.data[ 2]*(
b(0,2) +
b(2,0)) + a.data[ 3]*
b(1,1) + a.data[ 4]*(
b(1,2) +
b(2,1)) + a.data[ 5]*
b(2,2);
680 ret(1,1) += a.data[ 3]*
b(0,0) + a.data[ 8]*(
b(0,1) +
b(1,0)) + a.data[ 12]*(
b(0,2) +
b(2,0)) + a.data[ 15]*
b(1,1) + a.data[ 16]*(
b(1,2) +
b(2,1)) + a.data[ 17]*
b(2,2);
681 ret(2,2) += a.data[ 5]*
b(0,0) + a.data[ 10]*(
b(0,1) +
b(1,0)) + a.data[ 14]*(
b(0,2) +
b(2,0)) + a.data[ 17]*
b(1,1) + a.data[ 19]*(
b(1,2) +
b(2,1)) + a.data[ 20]*
b(2,2);
682 ret(1,2) += a.data[ 4]*
b(0,0) + a.data[ 9]*(
b(0,1) +
b(1,0)) + a.data[ 13]*(
b(0,2) +
b(2,0)) + a.data[ 16]*
b(1,1) + a.data[ 18]*(
b(1,2) +
b(2,1)) + a.data[ 19]*
b(2,2);
683 ret(2,0) += a.data[ 2]*
b(0,0) + a.data[ 7]*(
b(0,1) +
b(1,0)) + a.data[ 11]*(
b(0,2) +
b(2,0)) + a.data[ 12]*
b(1,1) + a.data[ 13]*(
b(1,2) +
b(2,1)) + a.data[ 14]*
b(2,2);
684 ret(0,1) += a.data[ 1]*
b(0,0) + a.data[ 6]*(
b(0,1) +
b(1,0)) + a.data[ 7]*(
b(0,2) +
b(2,0)) + a.data[ 8]*
b(1,1) + a.data[ 9]*(
b(1,2) +
b(2,1)) + a.data[ 10]*
b(2,2);
695 Eigen::Matrix<amrex::Real,2,2>
ret = Eigen::Matrix<amrex::Real,2,2>::Zero();
696 for (
int i = 0;
i < 2;
i++)
697 for (
int j = 0;
j < 2;
j++)
698 for (
int k = 0;
k < 2;
k++)
699 for (
int l = 0;
l < 2;
l++)
716 a.data[ 0]*
b(0,0,0) +
717 a.data[ 1]*(
b(0,1,0) +
b(1,0,0) +
b(0,0,1)) +
718 a.data[ 2]*(
b(0,2,0) +
b(2,0,0) +
b(0,0,2)) +
719 a.data[ 3]*
b(1,1,0) +
720 a.data[ 4]*(
b(1,2,0) +
b(2,1,0)) +
721 a.data[ 5]*
b(2,2,0) +
722 a.data[ 6]*(
b(0,1,1) +
b(1,0,1)) +
723 a.data[ 7]*(
b(0,2,1) +
b(2,0,1) +
b(0,1,2) +
b(1,0,2)) +
724 a.data[ 8]*
b(1,1,1) +
725 a.data[ 9]*(
b(1,2,1) +
b(2,1,1)) +
726 a.data[ 10]*
b(2,2,1) +
727 a.data[ 11]*(
b(0,2,2) +
b(2,0,2)) +
728 a.data[ 12]*
b(1,1,2) +
729 a.data[ 13]*(
b(1,2,2) +
b(2,1,2)) +
730 a.data[ 14]*
b(2,2,2);
733 a.data[ 1]*
b(0,0,0) +
734 a.data[ 4]*
b(0,0,2) +
735 a.data[ 3]*
b(0,0,1) +
736 a.data[ 6]*(
b(0,1,0) +
b(1,0,0)) +
737 a.data[ 7]*(
b(0,2,0) +
b(2,0,0)) +
738 a.data[ 8]*(
b(1,1,0) +
b(0,1,1) +
b(1,0,1)) +
739 a.data[ 9]*(
b(1,2,0) +
b(2,1,0) +
b(0,1,2) +
b(1,0,2)) +
740 a.data[ 10]*
b(2,2,0) +
741 a.data[ 12]*(
b(0,2,1) +
b(2,0,1)) +
742 a.data[ 13]*(
b(0,2,2) +
b(2,0,2)) +
743 a.data[ 15]*
b(1,1,1) +
744 a.data[ 16]*(
b(1,2,1) +
b(2,1,1) +
b(1,1,2)) +
745 a.data[ 17]*
b(2,2,1) +
746 a.data[ 18]*(
b(1,2,2) +
b(2,1,2)) +
747 a.data[ 19]*
b(2,2,2);
750 a.data[ 2]*
b(0,0,0) +
751 a.data[ 4]*
b(0,0,1) +
752 a.data[ 5]*
b(0,0,2) +
753 a.data[ 7]*(
b(0,1,0) +
b(1,0,0)) +
754 a.data[ 9]*(
b(0,1,1) +
b(1,0,1)) +
755 a.data[ 10]*(
b(0,1,2) +
b(1,0,2)) +
756 a.data[ 11]*(
b(0,2,0) +
b(2,0,0)) +
757 a.data[ 12]*
b(1,1,0) +
758 a.data[ 13]*(
b(1,2,0) +
b(2,1,0) +
b(0,2,1) +
b(2,0,1)) +
759 a.data[ 14]*(
b(2,2,0) +
b(0,2,2) +
b(2,0,2)) +
760 a.data[ 16]*
b(1,1,1) +
761 a.data[ 17]*
b(1,1,2) +
762 a.data[ 18]*(
b(1,2,1) +
b(2,1,1)) +
763 a.data[ 19]*(
b(2,2,1) +
b(1,2,2) +
b(2,1,2)) +
764 a.data[ 20]*
b(2,2,2);
Matrix< _Scalar, _Rows, _Cols > & operator*=(const amrex::Vector< amrex::Real > &x)
AMREX_FORCE_INLINE void operator+=(const OP_CLASS &rhs)
static Matrix4< 2, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 2, Sym::MajorMinor > Random()
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 2, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
static Matrix4< 2, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
static Matrix4< 2, Sym::MajorMinor > Zero()
static Matrix4< 2, Sym::MajorMinor > Increment()
void Print(std::ostream &os)
bool contains_nan() const
AMREX_GPU_HOST_DEVICE Matrix4()
static Matrix4< 2, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 3, Sym::MajorMinor > Increment()
static Matrix4< 3, Sym::MajorMinor > Random()
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 3, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
AMREX_GPU_HOST_DEVICE Matrix4()
void Print(std::ostream &os)
bool contains_nan() const
static Matrix4< 3, Sym::MajorMinor > Zero()
static Matrix4< 3, Sym::MajorMinor > Transverse(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
static Matrix4< 3, Sym::MajorMinor > Trans(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Quaternion q)
static Matrix4< 3, Sym::MajorMinor > Cubic(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R=Eigen::Matrix3d::Identity())
A collection of data types and symmetry-reduced data structures.
AMREX_FORCE_INLINE Quaternion operator-(const Quaternion a, const Quaternion b)
AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Matrix4< AMREX_SPACEDIM, Sym::Diagonal > operator/(const Matrix4< AMREX_SPACEDIM, Sym::Diagonal > &a, const Set::Scalar &b)
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
AMREX_FORCE_INLINE bool operator==(const Quaternion a, const Quaternion b)
AMREX_FORCE_INLINE Quaternion operator+(const Quaternion a, const Quaternion b)
AMREX_FORCE_INLINE Quaternion operator*(const Set::Scalar alpha, const Quaternion b)
void Abort(const char *msg)