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