1#ifndef SET_MATRIX4_MAJORMINOR_H
2#define SET_MATRIX4_MAJORMINOR_H
28 if (
uid==0 )
return data[0];
29 else if (
uid==8 ||
uid==4 ||
uid==2 ||
uid==1 )
return data[1];
30 else if (
uid==12 ||
uid==3 )
return data[2];
31 else if (
uid==10 ||
uid==6 ||
uid==9 ||
uid==5 )
return data[3];
32 else if (
uid==14 ||
uid==13 ||
uid==11 ||
uid==7 )
return data[4];
33 else if (
uid==15 )
return data[5];
41 if (
uid==0 )
return data[0];
42 else if (
uid==8 ||
uid==4 ||
uid==2 ||
uid==1 )
return data[1];
43 else if (
uid==12 ||
uid==3 )
return data[2];
44 else if (
uid==10 ||
uid==6 ||
uid==9 ||
uid==5 )
return data[3];
45 else if (
uid==14 ||
uid==13 ||
uid==11 ||
uid==7 )
return data[4];
46 else if (
uid==15 )
return data[5];
53 for (
int k = 0;
k < 2;
k++)
55 for (
int i = -1;
i < 3;
i++)
57 for (
int l = 0;
l < 2;
l++)
59 if (
i==-1)
os <<
"┌ ┐ ";
60 else if (
i == 2)
os <<
"└ ┘ ";
64 for (
int j = 0;
j < 2;
j++)
67 os << std::scientific << std::setw(11) << val ;
80 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]);
84 if (std::isnan(data[ 0]))
return true;
85 if (std::isnan(data[ 1]))
return true;
86 if (std::isnan(data[ 2]))
return true;
87 if (std::isnan(data[ 3]))
return true;
88 if (std::isnan(data[ 4]))
return true;
89 if (std::isnan(data[ 5]))
return true;
119 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
124 for(
int i = 0;
i < 3;
i++)
125 for(
int j = 0;
j < 3;
j++)
126 for(
int k = 0;
k < 3;
k++)
127 for(
int l = 0;
l < 3;
l++)
134 for(
int p = 0;
p < 2;
p++)
135 for(
int q = 0; q < 2; q++)
136 for(
int s = 0;
s < 2;
s++)
137 for(
int t = 0;
t < 2;
t++)
140 for(
int i = 0;
i < 3;
i++)
141 for(
int j = 0;
j < 3;
j++)
142 for(
int k = 0;
k < 3;
k++)
143 for(
int l = 0;
l < 3;
l++)
152 R = Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
153 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
154 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
160 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
161 return Cubic(C11,C12,C44,R);
170 for(
int i = 0;
i < 3;
i++)
171 for(
int j = 0;
j < 3;
j++)
172 for(
int k = 0;
k < 3;
k++)
173 for(
int l = 0;
l < 3;
l++)
180 else if(
i==
j &&
k==
l)
185 else if((
i==
k &&
j==
l) && (
i!=
j))
196 Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
197 Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
199 for(
int p = 0;
p < 2;
p++)
200 for(
int q = 0; q < 2; q++)
201 for(
int s = 0;
s < 2;
s++)
202 for(
int t = 0;
t < 2;
t++)
205 for(
int i = 0;
i < 3;
i++)
206 for(
int j = 0;
j < 3;
j++)
207 for(
int k = 0;
k < 3;
k++)
208 for(
int l = 0;
l < 3;
l++)
215 Eigen::Quaterniond q =
216 Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
217 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
218 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
219 Eigen::Matrix3d R = q.toRotationMatrix();
224 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
225 return Transverse(C11,C12,
C13,
C33,C44,R);
241 for (
int i = 0 ;
i < 6 ;
i++)
if (a.data[
i] !=
b.data[
i])
return false;
248 for (
int i = 0 ;
i < 6 ;
i++)
ret.data[
i] = a.data[
i] +
b.data[
i];
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;
274 for (
int i = 0 ;
i < 6 ;
i++)
ret.data[
i] = a.data[
i] /
b;
280 Eigen::Matrix<Set::Scalar,2,2>
ret = Eigen::Matrix<Set::Scalar,2,2>::Zero();
281 ret(0,0) = a.data[0]*
b(0,0) + a.data[1]*(
b(0,1) +
b(1,0)) + a.data[2]*
b(1,1);
282 ret(0,1) = a.data[1]*
b(0,0) + a.data[3]*(
b(0,1) +
b(1,0)) + a.data[4]*
b(1,1);
283 ret(1,1) = a.data[2]*
b(0,0) + a.data[4]*(
b(0,1) +
b(1,0)) + a.data[5]*
b(1,1);
294 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);
295 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);
323 int uid =
i + 3*
j + 9*
k + 27*
l;
325 if (
uid==0 )
return data[0];
327 else if (
uid==27 ||
uid==9 ||
uid==3 ||
uid==1 )
return data[1];
329 else if (
uid==54 ||
uid==18 ||
uid==6 ||
uid==2 )
return data[2];
331 else if (
uid==36 ||
uid==4 )
return data[3];
333 else if (
uid==63 ||
uid==45 ||
uid==7 ||
uid==5 )
return data[4];
335 else if (
uid==72 ||
uid==8 )
return data[5];
337 else if (
uid==30 ||
uid==12 ||
uid==28 ||
uid==10 )
return data[6];
341 else if (
uid==39 ||
uid==37 ||
uid==31 ||
uid==13 )
return data[8];
345 else if (
uid==75 ||
uid==73 ||
uid==35 ||
uid==17 )
return data[10];
347 else if (
uid==60 ||
uid==24 ||
uid==56 ||
uid==20 )
return data[11];
349 else if (
uid==42 ||
uid==58 ||
uid==22 ||
uid==38 )
return data[12];
353 else if (
uid==78 ||
uid==74 ||
uid==62 ||
uid==26 )
return data[14];
355 else if (
uid==40 )
return data[15];
357 else if (
uid==67 ||
uid==49 ||
uid==43 ||
uid==41 )
return data[16];
359 else if (
uid==76 ||
uid==44 )
return data[17];
361 else if (
uid==70 ||
uid==52 ||
uid==68 ||
uid==50 )
return data[18];
363 else if (
uid==79 ||
uid==77 ||
uid==71 ||
uid==53 )
return data[19];
365 else if (
uid==80 )
return data[20];
372 int uid =
i + 3*
j + 9*
k + 27*
l;
374 if (
uid==0 )
return data[0];
376 else if (
uid==27 ||
uid==9 ||
uid==3 ||
uid==1 )
return data[1];
378 else if (
uid==54 ||
uid==18 ||
uid==6 ||
uid==2 )
return data[2];
380 else if (
uid==36 ||
uid==4 )
return data[3];
382 else if (
uid==63 ||
uid==45 ||
uid==7 ||
uid==5 )
return data[4];
384 else if (
uid==72 ||
uid==8 )
return data[5];
386 else if (
uid==30 ||
uid==12 ||
uid==28 ||
uid==10 )
return data[6];
390 else if (
uid==39 ||
uid==37 ||
uid==31 ||
uid==13 )
return data[8];
394 else if (
uid==75 ||
uid==73 ||
uid==35 ||
uid==17 )
return data[10];
396 else if (
uid==60 ||
uid==24 ||
uid==56 ||
uid==20 )
return data[11];
398 else if (
uid==42 ||
uid==58 ||
uid==22 ||
uid==38 )
return data[12];
402 else if (
uid==78 ||
uid==74 ||
uid==62 ||
uid==26 )
return data[14];
404 else if (
uid==40 )
return data[15];
406 else if (
uid==67 ||
uid==49 ||
uid==43 ||
uid==41 )
return data[16];
408 else if (
uid==76 ||
uid==44 )
return data[17];
410 else if (
uid==70 ||
uid==52 ||
uid==68 ||
uid==50 )
return data[18];
412 else if (
uid==79 ||
uid==77 ||
uid==71 ||
uid==53 )
return data[19];
414 else if (
uid==80 )
return data[20];
420 for (
int i = 0;
i < 21;
i++)
421 os <<
"i = " <<
i <<
" " << data[
i] << std::endl;
424 for (
int k = 0;
k < 3;
k++)
426 for (
int i = -1;
i < 4;
i++)
428 for (
int l = 0;
l < 3;
l++)
430 if (
i==-1)
os <<
"┌ ┐ ";
431 else if (
i == 3)
os <<
"└ ┘ ";
435 for (
int j = 0;
j < 3;
j++)
438 os << std::scientific << std::setw(11) << val ;
452 for (
int i = 0;
i < 21;
i++)
retsq += data[
i];
453 return std::sqrt(
retsq);
457 if (std::isnan(data[ 0]))
return true;
458 if (std::isnan(data[ 1]))
return true;
459 if (std::isnan(data[ 2]))
return true;
460 if (std::isnan(data[ 3]))
return true;
461 if (std::isnan(data[ 4]))
return true;
462 if (std::isnan(data[ 5]))
return true;
463 if (std::isnan(data[ 6]))
return true;
464 if (std::isnan(data[ 7]))
return true;
465 if (std::isnan(data[ 8]))
return true;
466 if (std::isnan(data[ 9]))
return true;
467 if (std::isnan(data[10]))
return true;
468 if (std::isnan(data[11]))
return true;
469 if (std::isnan(data[12]))
return true;
470 if (std::isnan(data[13]))
return true;
471 if (std::isnan(data[14]))
return true;
472 if (std::isnan(data[15]))
return true;
473 if (std::isnan(data[16]))
return true;
474 if (std::isnan(data[17]))
return true;
475 if (std::isnan(data[18]))
return true;
476 if (std::isnan(data[19]))
return true;
477 if (std::isnan(data[20]))
return true;
508 Eigen::Matrix3d R = Eigen::Matrix3d::Identity())
513 for(
int i = 0;
i < 3;
i++)
514 for(
int j = 0;
j < 3;
j++)
515 for(
int k = 0;
k < 3;
k++)
516 for(
int l = 0;
l < 3;
l++)
523 for(
int p = 0;
p < 3;
p++)
524 for(
int q = 0; q < 3; q++)
525 for(
int s = 0;
s < 3;
s++)
526 for(
int t = 0;
t < 3;
t++)
529 for(
int i = 0;
i < 3;
i++)
530 for(
int j = 0;
j < 3;
j++)
531 for(
int k = 0;
k < 3;
k++)
532 for(
int l = 0;
l < 3;
l++)
541 R = Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
542 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
543 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
549 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
550 return Cubic(C11,C12,C44,R);
558 for(
int i = 0;
i < 3;
i++)
559 for(
int j = 0;
j < 3;
j++)
560 for(
int k = 0;
k < 3;
k++)
561 for(
int l = 0;
l < 3;
l++)
568 else if(
i==
j &&
k==
l)
573 else if((
i==
k &&
j==
l) && (
i!=
j))
584 Ctmp[0][1][0][1] = 0.5 * (C11 - C12);
585 Ctmp[1][0][1][0] = 0.5 * (C11 - C12);
587 for(
int p = 0;
p < 3;
p++)
588 for(
int q = 0; q < 3; q++)
589 for(
int s = 0;
s < 3;
s++)
590 for(
int t = 0;
t < 3;
t++)
593 for(
int i = 0;
i < 3;
i++)
594 for(
int j = 0;
j < 3;
j++)
595 for(
int k = 0;
k < 3;
k++)
596 for(
int l = 0;
l < 3;
l++)
603 Eigen::Quaterniond q =
604 Eigen::AngleAxisd(
phi2, Eigen::Vector3d::UnitX()) *
605 Eigen::AngleAxisd(
Phi, Eigen::Vector3d::UnitZ()) *
606 Eigen::AngleAxisd(
phi1, Eigen::Vector3d::UnitX());
607 Eigen::Matrix3d R = q.toRotationMatrix();
612 Eigen::Matrix3d R = q.normalized().toRotationMatrix();
613 return Transverse(C11,C12,
C13,
C33,C44,R);
628 for (
int i = 0 ;
i < 21 ;
i++)
if (a.data[
i] !=
b.data[
i])
return false;
635 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] +
b.data[
i];
642 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] -
b.data[
i];
649 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] *
b;
661 for (
int i = 0 ;
i < 21 ;
i++)
ret.data[
i] = a.data[
i] /
b;
667 Eigen::Matrix<amrex::Real,3,3>
ret = Eigen::Matrix<amrex::Real,3,3>::Zero();
668 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);
669 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);
670 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);
671 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);
672 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);
673 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);
684 Eigen::Matrix<amrex::Real,2,2>
ret = Eigen::Matrix<amrex::Real,2,2>::Zero();
685 for (
int i = 0;
i < 2;
i++)
686 for (
int j = 0;
j < 2;
j++)
687 for (
int k = 0;
k < 2;
k++)
688 for (
int l = 0;
l < 2;
l++)
705 a.data[ 0]*
b(0,0,0) +
706 a.data[ 1]*(
b(0,1,0) +
b(1,0,0) +
b(0,0,1)) +
707 a.data[ 2]*(
b(0,2,0) +
b(2,0,0) +
b(0,0,2)) +
708 a.data[ 3]*
b(1,1,0) +
709 a.data[ 4]*(
b(1,2,0) +
b(2,1,0)) +
710 a.data[ 5]*
b(2,2,0) +
711 a.data[ 6]*(
b(0,1,1) +
b(1,0,1)) +
712 a.data[ 7]*(
b(0,2,1) +
b(2,0,1) +
b(0,1,2) +
b(1,0,2)) +
713 a.data[ 8]*
b(1,1,1) +
714 a.data[ 9]*(
b(1,2,1) +
b(2,1,1)) +
715 a.data[ 10]*
b(2,2,1) +
716 a.data[ 11]*(
b(0,2,2) +
b(2,0,2)) +
717 a.data[ 12]*
b(1,1,2) +
718 a.data[ 13]*(
b(1,2,2) +
b(2,1,2)) +
719 a.data[ 14]*
b(2,2,2);
722 a.data[ 1]*
b(0,0,0) +
723 a.data[ 4]*
b(0,0,2) +
724 a.data[ 3]*
b(0,0,1) +
725 a.data[ 6]*(
b(0,1,0) +
b(1,0,0)) +
726 a.data[ 7]*(
b(0,2,0) +
b(2,0,0)) +
727 a.data[ 8]*(
b(1,1,0) +
b(0,1,1) +
b(1,0,1)) +
728 a.data[ 9]*(
b(1,2,0) +
b(2,1,0) +
b(0,1,2) +
b(1,0,2)) +
729 a.data[ 10]*
b(2,2,0) +
730 a.data[ 12]*(
b(0,2,1) +
b(2,0,1)) +
731 a.data[ 13]*(
b(0,2,2) +
b(2,0,2)) +
732 a.data[ 15]*
b(1,1,1) +
733 a.data[ 16]*(
b(1,2,1) +
b(2,1,1) +
b(1,1,2)) +
734 a.data[ 17]*
b(2,2,1) +
735 a.data[ 18]*(
b(1,2,2) +
b(2,1,2)) +
736 a.data[ 19]*
b(2,2,2);
739 a.data[ 2]*
b(0,0,0) +
740 a.data[ 4]*
b(0,0,1) +
741 a.data[ 5]*
b(0,0,2) +
742 a.data[ 7]*(
b(0,1,0) +
b(1,0,0)) +
743 a.data[ 9]*(
b(0,1,1) +
b(1,0,1)) +
744 a.data[ 10]*(
b(0,1,2) +
b(1,0,2)) +
745 a.data[ 11]*(
b(0,2,0) +
b(2,0,0)) +
746 a.data[ 12]*
b(1,1,0) +
747 a.data[ 13]*(
b(1,2,0) +
b(2,1,0) +
b(0,2,1) +
b(2,0,1)) +
748 a.data[ 14]*(
b(2,2,0) +
b(0,2,2) +
b(2,0,2)) +
749 a.data[ 16]*
b(1,1,1) +
750 a.data[ 17]*
b(1,1,2) +
751 a.data[ 18]*(
b(1,2,1) +
b(2,1,1)) +
752 a.data[ 19]*(
b(2,2,1) +
b(1,2,2) +
b(2,1,2)) +
753 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 > 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 > Randomize()
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 > Randomize()
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)