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 : #include "Util/Util.H"
16 :
17 : namespace Set
18 : {
19 : using Scalar = amrex::Real;
20 : using Vector = Eigen::Matrix<amrex::Real,AMREX_SPACEDIM,1>;
21 : using Vector3d = Eigen::Vector3d;
22 : using Covector = Eigen::Matrix<amrex::Real,1,AMREX_SPACEDIM>;
23 : using Matrix = Eigen::Matrix<amrex::Real,AMREX_SPACEDIM,AMREX_SPACEDIM>;
24 : using Matrix3d = Eigen::Matrix3d;
25 : using iMatrix = Eigen::Matrix<int,AMREX_SPACEDIM,AMREX_SPACEDIM>;
26 :
27 : AMREX_FORCE_INLINE
28 : Set::Matrix reduce(const Set::Matrix3d &in)
29 : {
30 1 : return in.block<AMREX_SPACEDIM,AMREX_SPACEDIM>(0,0);
31 : }
32 :
33 : AMREX_FORCE_INLINE
34 : Set::Matrix3d expand(const Set::Matrix &in)
35 : {
36 0 : Set::Matrix3d ret = Set::Matrix3d::Zero();
37 0 : ret.block<AMREX_SPACEDIM,AMREX_SPACEDIM>(0,0) = in;
38 0 : return ret;
39 : }
40 :
41 : //using Quaternion = Eigen::Quaterniond;
42 :
43 : class Quaternion : public Eigen::Quaterniond
44 : {
45 : public:
46 8 : Quaternion()
47 8 : : Eigen::Quaterniond() {}
48 144 : Quaternion(Set::Scalar w, Set::Scalar x, Set::Scalar y, Set::Scalar z)
49 144 : : Eigen::Quaterniond(w,x,y,z) {}
50 0 : Quaternion(const Eigen::Matrix3d &R)
51 0 : : Eigen::Quaterniond(R) {}
52 : Quaternion(Eigen::Quaterniond &q)
53 : : Eigen::Quaterniond(q.w(),q.x(),q.y(),q.z()) {}
54 :
55 : AMREX_FORCE_INLINE
56 : void operator = (const Eigen::Quaterniond &rhs)
57 : {
58 128 : w() = rhs.w(); x() = rhs.x(); y() = rhs.y(); z() = rhs.z();
59 128 : }
60 :
61 : AMREX_FORCE_INLINE
62 : void operator += (const Quaternion &rhs)
63 : {
64 0 : w() += rhs.w(); x() += rhs.x(); y() += rhs.y(); z() += rhs.z();
65 0 : }
66 : // AMREX_FORCE_INLINE
67 : // void operator * (const Set::Scalar alpha)
68 : // {
69 : // w += alpha; x *= alpha; y *= alpha; z *= alpha;
70 : // }
71 : friend Quaternion operator * (const Set::Scalar alpha, const Quaternion b);
72 : friend Quaternion operator * (const Quaternion b,const Set::Scalar alpha);
73 : friend Quaternion operator + (const Quaternion a, const Quaternion b);
74 : friend Quaternion operator - (const Quaternion a, const Quaternion b);
75 : friend bool operator == (const Quaternion a, const Quaternion b);
76 : };
77 : AMREX_FORCE_INLINE
78 : Quaternion operator * (const Set::Scalar alpha, const Quaternion b)
79 : {
80 : Quaternion ret;
81 : ret.w() = b.w()*alpha; ret.x() = b.x()*alpha; ret.y() = b.y()*alpha; ret.z() = b.z()*alpha;
82 : return ret;
83 : }
84 : AMREX_FORCE_INLINE
85 : Quaternion operator * (const Quaternion b, const Set::Scalar alpha)
86 : {
87 8 : Quaternion ret;
88 8 : ret.w() = b.w()*alpha; ret.x() = b.x()*alpha; ret.y() = b.y()*alpha; ret.z() = b.z()*alpha;
89 8 : return ret;
90 : }
91 : AMREX_FORCE_INLINE
92 : Quaternion operator + (const Quaternion a, const Quaternion b)
93 : {
94 0 : Quaternion ret;
95 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();
96 0 : return ret;
97 : }
98 :
99 : AMREX_FORCE_INLINE
100 : Quaternion operator - (const Quaternion a, const Quaternion b)
101 : {
102 0 : Quaternion ret;
103 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();
104 0 : return ret;
105 : }
106 :
107 : AMREX_FORCE_INLINE
108 : bool operator == (const Quaternion a, const Quaternion b)
109 : {
110 20 : if (fabs((a.w() - b.w())/(a.w()+a.w())) > 1E-8) return false;
111 20 : if (fabs((a.x() - b.x())/(a.x()+a.x())) > 1E-8) return false;
112 20 : if (fabs((a.y() - b.y())/(a.y()+a.y())) > 1E-8) return false;
113 20 : if (fabs((a.z() - b.z())/(a.z()+a.z())) > 1E-8) return false;
114 20 : return true;
115 : }
116 :
117 :
118 : static Scalar Garbage = NAN;
119 :
120 : AMREX_FORCE_INLINE
121 : Vector Position(const int &i, const int &j, const int &k, const amrex::Geometry &geom, const amrex::IndexType &ixType)
122 : {
123 : (void)k;
124 748028 : if (ixType == amrex::IndexType::TheNodeType())
125 : {
126 147740 : return Vector(AMREX_D_DECL(
127 : geom.ProbLo()[0] + ((amrex::Real)(i)) * geom.CellSize()[0],
128 : geom.ProbLo()[1] + ((amrex::Real)(j)) * geom.CellSize()[1],
129 73870 : geom.ProbLo()[2] + ((amrex::Real)(k)) * geom.CellSize()[2]));
130 : }
131 600288 : else if (ixType == amrex::IndexType::TheCellType())
132 : {
133 600288 : return Vector(AMREX_D_DECL(
134 : geom.ProbLo()[0] + ((amrex::Real)(i) + 0.5) * geom.CellSize()[0],
135 : geom.ProbLo()[1] + ((amrex::Real)(j) + 0.5) * geom.CellSize()[1],
136 300144 : geom.ProbLo()[2] + ((amrex::Real)(k) + 0.5) * geom.CellSize()[2]));
137 : }
138 0 : Util::Abort(INFO,"Unsupported index type");
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 : Util::Abort(INFO,"Unsupported index type");
194 : return Set::Vector::Zero();
195 : }
196 :
197 : enum Sym {None,Major,Minor,MajorMinor,Diagonal,Full,Isotropic};
198 : template<int dim,int sym> class Matrix4{};
199 : template<int dim, int sym>
200 : std::ostream&
201 0 : operator<< (std::ostream& os, const Matrix4<dim,sym>& b)
202 : {
203 0 : Matrix4<dim,sym> bcopy = b;
204 0 : bcopy.Print(os);
205 0 : return os;
206 : }
207 : }
208 : #endif
|
