Line data Source code
1 : #ifndef SET_BASE_H
2 : #define SET_BASE_H
3 :
4 : #include "AMReX.H"
5 : #include "AMReX_REAL.H"
6 : #include "AMReX_SPACE.H"
7 : #include "AMReX_Vector.H"
8 : #include "AMReX_BLassert.H"
9 : #include "AMReX_Geometry.H"
10 : #include "AMReX_IndexType.H"
11 : #include <eigen3/Eigen/Core>
12 : #include <eigen3/Eigen/Geometry>
13 : #include <eigen3/Eigen/SVD>
14 : #include <AMReX_REAL.H>
15 :
16 : namespace Set
17 : {
18 : using Scalar = amrex::Real;
19 : using Vector = Eigen::Matrix<amrex::Real,AMREX_SPACEDIM,1>;
20 : using Vector3d = Eigen::Vector3d;
21 : using Covector = Eigen::Matrix<amrex::Real,1,AMREX_SPACEDIM>;
22 : using Matrix = Eigen::Matrix<amrex::Real,AMREX_SPACEDIM,AMREX_SPACEDIM>;
23 : using Matrix3d = Eigen::Matrix3d;
24 : using iMatrix = Eigen::Matrix<int,AMREX_SPACEDIM,AMREX_SPACEDIM>;
25 :
26 : AMREX_FORCE_INLINE
27 : Set::Matrix reduce(const Set::Matrix3d &in)
28 : {
29 1 : return in.block<AMREX_SPACEDIM,AMREX_SPACEDIM>(0,0);
30 : }
31 :
32 : AMREX_FORCE_INLINE
33 : Set::Matrix3d expand(const Set::Matrix &in)
34 : {
35 0 : Set::Matrix3d ret = Set::Matrix3d::Zero();
36 0 : ret.block<AMREX_SPACEDIM,AMREX_SPACEDIM>(0,0) = in;
37 0 : return ret;
38 : }
39 :
40 : //using Quaternion = Eigen::Quaterniond;
41 :
42 : class Quaternion : public Eigen::Quaterniond
43 : {
44 : public:
45 8 : Quaternion()
46 8 : : Eigen::Quaterniond() {}
47 144 : Quaternion(Set::Scalar w, Set::Scalar x, Set::Scalar y, Set::Scalar z)
48 144 : : Eigen::Quaterniond(w,x,y,z) {}
49 0 : Quaternion(const Eigen::Matrix3d &R)
50 0 : : Eigen::Quaterniond(R) {}
51 : Quaternion(Eigen::Quaterniond &q)
52 : : Eigen::Quaterniond(q.w(),q.x(),q.y(),q.z()) {}
53 :
54 : AMREX_FORCE_INLINE
55 : void operator = (const Eigen::Quaterniond &rhs)
56 : {
57 128 : w() = rhs.w(); x() = rhs.x(); y() = rhs.y(); z() = rhs.z();
58 128 : }
59 :
60 : AMREX_FORCE_INLINE
61 : void operator += (const Quaternion &rhs)
62 : {
63 0 : w() += rhs.w(); x() += rhs.x(); y() += rhs.y(); z() += rhs.z();
64 0 : }
65 : // AMREX_FORCE_INLINE
66 : // void operator * (const Set::Scalar alpha)
67 : // {
68 : // w += alpha; x *= alpha; y *= alpha; z *= alpha;
69 : // }
70 : friend Quaternion operator * (const Set::Scalar alpha, const Quaternion b);
71 : friend Quaternion operator * (const Quaternion b,const Set::Scalar alpha);
72 : friend Quaternion operator + (const Quaternion a, const Quaternion b);
73 : friend Quaternion operator - (const Quaternion a, const Quaternion b);
74 : friend bool operator == (const Quaternion a, const Quaternion b);
75 : };
76 : AMREX_FORCE_INLINE
77 : Quaternion operator * (const Set::Scalar alpha, const Quaternion b)
78 : {
79 : Quaternion ret;
80 : ret.w() = b.w()*alpha; ret.x() = b.x()*alpha; ret.y() = b.y()*alpha; ret.z() = b.z()*alpha;
81 : return ret;
82 : }
83 : AMREX_FORCE_INLINE
84 : Quaternion operator * (const Quaternion b, const Set::Scalar alpha)
85 : {
86 8 : Quaternion ret;
87 8 : ret.w() = b.w()*alpha; ret.x() = b.x()*alpha; ret.y() = b.y()*alpha; ret.z() = b.z()*alpha;
88 8 : return ret;
89 : }
90 : AMREX_FORCE_INLINE
91 : Quaternion operator + (const Quaternion a, const Quaternion b)
92 : {
93 0 : Quaternion ret;
94 0 : ret.w() = a.w()+b.w(); ret.x() = a.x()+b.x(); ret.y() = a.y()+b.y(); ret.z() = a.z()+b.z();
95 0 : return ret;
96 : }
97 :
98 : AMREX_FORCE_INLINE
99 : Quaternion operator - (const Quaternion a, const Quaternion b)
100 : {
101 0 : Quaternion ret;
102 0 : ret.w() = a.w()-b.w(); ret.x() = a.x()-b.x(); ret.y() = a.y()-b.y(); ret.z() = a.z()-b.z();
103 0 : return ret;
104 : }
105 :
106 : AMREX_FORCE_INLINE
107 : bool operator == (const Quaternion a, const Quaternion b)
108 : {
109 20 : if (fabs((a.w() - b.w())/(a.w()+a.w())) > 1E-8) return false;
110 20 : if (fabs((a.x() - b.x())/(a.x()+a.x())) > 1E-8) return false;
111 20 : if (fabs((a.y() - b.y())/(a.y()+a.y())) > 1E-8) return false;
112 20 : if (fabs((a.z() - b.z())/(a.z()+a.z())) > 1E-8) return false;
113 20 : return true;
114 : }
115 :
116 :
117 : [[maybe_unused]] static Scalar Garbage = NAN;
118 :
119 : AMREX_FORCE_INLINE
120 : Vector Position(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
121 : {
122 : (void)k;
123 773798 : if (ixType == amrex::IndexType::TheNodeType())
124 : {
125 222710 : return Vector(AMREX_D_DECL(
126 : geom.ProbLo()[0] + ((amrex::Real)(i)) * geom.CellSize()[0],
127 : geom.ProbLo()[1] + ((amrex::Real)(j)) * geom.CellSize()[1],
128 111355 : geom.ProbLo()[2] + ((amrex::Real)(k)) * geom.CellSize()[2]));
129 : }
130 551088 : else if (ixType == amrex::IndexType::TheCellType())
131 : {
132 551088 : return Vector(AMREX_D_DECL(
133 : geom.ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom.CellSize()[0],
134 : geom.ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom.CellSize()[1],
135 275544 : geom.ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom.CellSize()[2]));
136 : }
137 0 : return Set::Vector::Zero();
138 : }
139 :
140 : AMREX_FORCE_INLINE
141 : Vector Normal( AMREX_D_DECL(bool xmin,bool ymin, bool zmin),
142 : AMREX_D_DECL(bool xmax,bool ymax, bool zmax))
143 : {
144 0 : Set::Vector N = Set::Vector::Zero();
145 0 : if (xmin) N(0) = -1.0;
146 0 : else if (xmax) N(0) = 1.0;
147 : #if AMREX_SPACEDIM>1
148 0 : if (ymin) N(1) = -1.0;
149 0 : else if (ymax) N(1) = 1.0;
150 : #endif
151 : #if AMREX_SPACEDIM>2
152 0 : if (zmin) N(2) = -1.0;
153 0 : else if (zmax) N(2) = 1.0;
154 : #endif
155 0 : return N;
156 : }
157 :
158 : AMREX_FORCE_INLINE
159 : Vector Size(const amrex::Geometry &geom)
160 : {
161 0 : Set::Vector size;
162 0 : size(0) = geom.ProbHi()[0] - geom.ProbLo()[0];
163 : #if AMREX_SPACEDIM>1
164 0 : size(1) = geom.ProbHi()[1] - geom.ProbLo()[1];
165 : #endif
166 : #if AMREX_SPACEDIM>2
167 0 : size(2) = geom.ProbHi()[2] - geom.ProbLo()[2];
168 : #endif
169 0 : return size;
170 : }
171 :
172 :
173 : AMREX_FORCE_INLINE
174 : Vector Volume(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
175 : {
176 : (void)k;
177 : if (ixType == amrex::IndexType::TheNodeType())
178 : {
179 : return Vector(AMREX_D_DECL(
180 : geom.ProbLo()[0] + ((amrex::Real)(i)) * geom.CellSize()[0],
181 : geom.ProbLo()[1] + ((amrex::Real)(j)) * geom.CellSize()[1],
182 : geom.ProbLo()[2] + ((amrex::Real)(k)) * geom.CellSize()[2]));
183 : }
184 : else if (ixType == amrex::IndexType::TheCellType())
185 : {
186 : return Vector(AMREX_D_DECL(
187 : geom.ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom.CellSize()[0],
188 : geom.ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom.CellSize()[1],
189 : geom.ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom.CellSize()[2]));
190 : }
191 : return Set::Vector::Zero();
192 : }
193 :
194 : enum Sym {None,Major,Minor,MajorMinor,Diagonal,Full,Isotropic};
195 : template<int dim,int sym> class Matrix4{};
196 : template<int dim, int sym>
197 : std::ostream&
198 0 : operator<< (std::ostream& os, const Matrix4<dim,sym>& b)
199 : {
200 0 : Matrix4<dim,sym> bcopy = b;
201 0 : bcopy.Print(os);
202 0 : return os;
203 : }
204 : }
205 : #endif
|
