36 #define LOKI_CLASS_LEVEL_THREADING
40 #define EPSILON 0.000001
42 #define CROSS(dest, v1, v2){ \
43 dest[0] = v1[1] * v2[2] - v1[2] * v2[1]; \
44 dest[1] = v1[2] * v2[0] - v1[0] * v2[2]; \
45 dest[2] = v1[0] * v2[1] - v1[1] * v2[0];}
47 #define DOT(v1, v2) (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2])
49 #define SUB(dest, v1, v2){ \
50 dest[0] = v1[0] - v2[0]; \
51 dest[1] = v1[1] - v2[1]; \
52 dest[2] = v1[2] - v2[2];}
68 Loki::SingleThreaded> Type;
69 return Type::Instance();
79 RiemannSolver::RiemannSolver() : m_requiresRotation(false),
116 int nFields = Fwd .num_elements();
117 int nPts = Fwd[0].num_elements();
122 for (
int i = 0; i < 3; ++i)
126 for (
int j = 0; j < nFields; ++j)
170 for (
int i = 0; i < inarray.num_elements(); ++i)
176 for (
int i = 0; i < vecLocs.num_elements(); i++)
178 ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
179 "vecLocs[i] element count mismatch");
181 switch (normals.num_elements())
189 const int nq = inarray[0].num_elements();
190 const int vx = (int)vecLocs[i][0];
191 const int vy = (int)vecLocs[i][1];
196 outarray[vx], 1, outarray[vx], 1);
200 outarray[vy], 1, outarray[vy], 1);
206 const int nq = inarray[0].num_elements();
207 const int vx = (int)vecLocs[i][0];
208 const int vy = (int)vecLocs[i][1];
209 const int vz = (int)vecLocs[i][2];
222 outarray[vx], 1, outarray[vx], 1);
227 outarray[vy], 1, outarray[vy], 1);
232 outarray[vz], 1, outarray[vz], 1);
237 ASSERTL1(
false,
"Invalid space dimension.");
257 for (
int i = 0; i < inarray.num_elements(); ++i)
263 for (
int i = 0; i < vecLocs.num_elements(); i++)
265 ASSERTL1(vecLocs[i].num_elements() == normals.num_elements(),
266 "vecLocs[i] element count mismatch");
268 switch (normals.num_elements())
276 const int nq = normals[0].num_elements();
277 const int vx = (int)vecLocs[i][0];
278 const int vy = (int)vecLocs[i][1];
283 outarray[vx], 1, outarray[vx], 1);
287 outarray[vy], 1, outarray[vy], 1);
293 const int nq = normals[0].num_elements();
294 const int vx = (int)vecLocs[i][0];
295 const int vy = (int)vecLocs[i][1];
296 const int vz = (int)vecLocs[i][2];
302 outarray[vx], 1, outarray[vx], 1);
307 outarray[vy], 1, outarray[vy], 1);
312 outarray[vz], 1, outarray[vz], 1);
317 ASSERTL1(
false,
"Invalid space dimension.");
397 const int nq = normals[0].num_elements();
404 for (i = 0; i < 9; ++i)
408 for (i = 0; i < normals[0].num_elements(); ++i)
412 tn[0] = normals[0][i];
413 tn[1] = normals[1][i];
414 tn[2] = normals[2][i];
417 for (j = 0; j < 9; ++j)
454 x[0] = (from[0] > 0.0)? from[0] : -from[0];
455 x[1] = (from[1] > 0.0)? from[1] : -from[1];
456 x[2] = (from[2] > 0.0)? from[2] : -from[2];
462 x[0] = 1.0; x[1] = x[2] = 0.0;
466 x[2] = 1.0; x[0] = x[1] = 0.0;
473 x[1] = 1.0; x[0] = x[2] = 0.0;
477 x[2] = 1.0; x[0] = x[1] = 0.0;
481 u[0] = x[0] - from[0];
482 u[1] = x[1] - from[1];
483 u[2] = x[2] - from[2];
484 v[0] = x[0] - to [0];
485 v[1] = x[1] - to [1];
486 v[2] = x[2] - to [2];
488 c1 = 2.0 /
DOT(u, u);
489 c2 = 2.0 /
DOT(v, v);
490 c3 = c1 * c2 *
DOT(u, v);
492 for (i = 0; i < 3; i++) {
493 for (j = 0; j < 3; j++) {
494 mat[3*i+j] = - c1 * u[i] * u[j]
510 mat[0] = e + hvx * v[0];
511 mat[1] = hvxy - v[2];
512 mat[2] = hvxz + v[1];
513 mat[3] = hvxy + v[2];
514 mat[4] = e + h * v[1] * v[1];
515 mat[5] = hvyz - v[0];
516 mat[6] = hvxz - v[1];
517 mat[7] = hvyz + v[0];
518 mat[8] = e + hvz * v[2];
void FromToRotation(Array< OneD, const NekDouble > &from, Array< OneD, const NekDouble > &to, NekDouble *mat)
A function for creating a rotation matrix that rotates a vector from into another vector to...
SOLVER_UTILS_EXPORT bool CheckScalars(std::string name)
Determine whether a scalar has been defined in m_scalars.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
SOLVER_UTILS_EXPORT bool CheckVectors(std::string name)
Determine whether a vector has been defined in m_vectors.
SOLVER_UTILS_EXPORT void rotateFromNormal(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
Rotate a vector field from trace normal.
virtual void v_Solve(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)=0
LibUtilities::NekFactory< std::string, RiemannSolver > RiemannSolverFactory
Datatype of the NekFactory used to instantiate classes derived from the RiemannSolver class...
#define CROSS(dest, v1, v2)
std::map< std::string, RSScalarFuncType > m_scalars
Map of scalar function types.
void GenerateRotationMatrices(const Array< OneD, const Array< OneD, NekDouble > > &normals)
Generate rotation matrices for 3D expansions.
RiemannSolverFactory & GetRiemannSolverFactory()
SOLVER_UTILS_EXPORT void Solve(const int nDim, const Array< OneD, const Array< OneD, NekDouble > > &Fwd, const Array< OneD, const Array< OneD, NekDouble > > &Bwd, Array< OneD, Array< OneD, NekDouble > > &flux)
Perform the Riemann solve given the forwards and backwards spaces.
std::map< std::string, RSScalarFuncType > m_auxScal
Map of auxiliary scalar function types.
Array< OneD, Array< OneD, Array< OneD, NekDouble > > > m_rotStorage
Rotation storage.
SOLVER_UTILS_EXPORT bool CheckAuxScal(std::string name)
Determine whether a scalar has been defined in m_auxScal.
Array< OneD, Array< OneD, NekDouble > > m_rotMat
Rotation matrices for each trace quadrature point.
SOLVER_UTILS_EXPORT bool CheckAuxVec(std::string name)
Determine whether a vector has been defined in m_auxVec.
std::map< std::string, RSVecFuncType > m_auxVec
Map of auxiliary vector function types.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
bool m_requiresRotation
Indicates whether the Riemann solver requires a rotation to be applied to the velocity fields...
void Vvtvm(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvm (vector times vector plus vector): z = w*x - y
std::map< std::string, RSParamFuncType > m_params
Map of parameter function types.
SOLVER_UTILS_EXPORT bool CheckParams(std::string name)
Determine whether a parameter has been defined in m_params.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
std::map< std::string, RSVecFuncType > m_vectors
Map of vector function types.
SOLVER_UTILS_EXPORT void rotateToNormal(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const Array< OneD, const Array< OneD, NekDouble > > &normals, const Array< OneD, const Array< OneD, NekDouble > > &vecLocs, Array< OneD, Array< OneD, NekDouble > > &outarray)
Rotate a vector field to trace normal.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Provides a generic Factory class.