Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in / Register
Toggle navigation
Menu
Open sidebar
Nektar
Nektar
Commits
4b560009
Commit
4b560009
authored
Apr 13, 2017
by
Spencer Sherwin
Browse files
Merge branch 'master' into fix/GetOffsetElmtId
parents
e5fa5e00
9738e4ae
Changes
97
Expand all
Hide whitespace changes
Inline
Side-by-side
CHANGELOG.md
View file @
4b560009
...
...
@@ -4,10 +4,13 @@ Changelog
v4.5.0
------
**NekMesh**
:
-
Add periodic boundary condition meshing in 2D (!733)
-
Adjust boundary layer thickness in corners in 2D (!739)
**Library**
-
Added in sum factorisation version for pyramid expansions and orthogonal
expansion in pyramids (!750)
**Documentation**
:
-
Added the developer-guide repository as a submodule (!751)
...
...
@@ -16,6 +19,20 @@ v4.4.1
**Library**
-
Remove m_offset_elmt_id and GetOffsetElmtId which fixed problems in 2D when
quad elements are listed before tri elements (!758)
-
Remove the duplicate output of errorutil (!756)
-
Fix interpolation issue with Lagrange basis functions (!768)
**FieldConvert**
:
-
Fix issue with FieldConvert when range flag used (!761)
**NekMesh**
:
-
Fix memory consumption issue with Gmsh output (!747, !762)
-
Rework meshing control so that if possible viewable meshes will be dumped
when some part of the system fails (!756)
-
Add manifold meshing option (!756)
**FieldConvert:**
-
Fix issue with field ordering in the interppointdatatofld module (!754)
v4.4.0
------
...
...
docs/doxygen/Doxyfile.in
View file @
4b560009
...
...
@@ -757,6 +757,7 @@ INPUT = @CMAKE_SOURCE_DIR@/docs/doxygen/ \
@CMAKE_SOURCE_DIR@/library/LibUtilities/ \
@CMAKE_SOURCE_DIR@/library/StdRegions/ \
@CMAKE_SOURCE_DIR@/library/SpatialDomains/ \
@CMAKE_SOURCE_DIR@/library/Collections/ \
@CMAKE_SOURCE_DIR@/library/LocalRegions/ \
@CMAKE_SOURCE_DIR@/library/MultiRegions/ \
@CMAKE_SOURCE_DIR@/library/GlobalMapping/ \
...
...
library/Collections/BwdTrans.cpp
View file @
4b560009
...
...
@@ -918,7 +918,7 @@ class BwdTrans_SumFac_Prism : public Operator
{
Blas
::
Dgemm
(
'N'
,
'N'
,
m_nquad2
,
m_numElmt
,
m_nmodes2
-
i
,
1.0
,
m_base2
.
get
()
+
mode
*
m_nquad2
,
m_nquad2
,
&
input
[
0
]
+
mode1
,
totmodes
,
0.0
,
input
.
get
()
+
mode1
,
totmodes
,
0.0
,
&
wsp
[
j
*
m_nquad2
*
m_numElmt
*
m_nmodes0
+
cnt
],
m_nquad2
);
mode1
+=
m_nmodes2
-
i
;
...
...
@@ -1018,5 +1018,168 @@ OperatorKey BwdTrans_SumFac_Prism::m_type = GetOperatorFactory().
OperatorKey
(
ePrism
,
eBwdTrans
,
eSumFac
,
false
),
BwdTrans_SumFac_Prism
::
create
,
"BwdTrans_SumFac_Prism"
);
/**
* @brief Backward transform operator using sum-factorisation (Pyr)
*/
class
BwdTrans_SumFac_Pyr
:
public
Operator
{
public:
OPERATOR_CREATE
(
BwdTrans_SumFac_Pyr
)
virtual
~
BwdTrans_SumFac_Pyr
()
{
}
virtual
void
operator
()(
const
Array
<
OneD
,
const
NekDouble
>
&
input
,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
output1
,
Array
<
OneD
,
NekDouble
>
&
output2
,
Array
<
OneD
,
NekDouble
>
&
wsp
)
{
ASSERTL1
(
wsp
.
num_elements
()
==
m_wspSize
,
"Incorrect workspace size"
);
// Assign second half of workspace for 2nd DGEMM operation.
int
totmodes
=
m_stdExp
->
GetNcoeffs
();
Array
<
OneD
,
NekDouble
>
wsp2
=
wsp
+
m_nmodes0
*
m_nmodes1
*
m_nquad2
*
m_numElmt
;
Vmath
::
Zero
(
m_nmodes0
*
m_nmodes1
*
m_nquad2
*
m_numElmt
,
wsp
,
1
);
int
i
=
0
;
int
j
=
0
;
int
mode
=
0
;
int
mode1
=
0
;
int
cnt
=
0
;
for
(
i
=
0
;
i
<
m_nmodes0
;
++
i
)
{
for
(
j
=
0
;
j
<
m_nmodes1
;
++
j
,
++
cnt
)
{
int
ijmax
=
max
(
i
,
j
);
Blas
::
Dgemm
(
'N'
,
'N'
,
m_nquad2
,
m_numElmt
,
m_nmodes2
-
ijmax
,
1.0
,
m_base2
.
get
()
+
mode
*
m_nquad2
,
m_nquad2
,
input
.
get
()
+
mode1
,
totmodes
,
0.0
,
wsp
.
get
()
+
cnt
*
m_nquad2
*
m_numElmt
,
m_nquad2
);
mode
+=
m_nmodes2
-
ijmax
;
mode1
+=
m_nmodes2
-
ijmax
;
}
//increment mode in case order1!=order2
for
(
j
=
m_nmodes1
;
j
<
m_nmodes2
-
i
;
++
j
)
{
int
ijmax
=
max
(
i
,
j
);
mode
+=
m_nmodes2
-
ijmax
;
}
}
// vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
// component is evaluated
if
(
m_sortTopVertex
)
{
for
(
i
=
0
;
i
<
m_numElmt
;
++
i
)
{
// top singular vertex
// (1+c)/2 x (1+b)/2 x (1-a)/2 component
Blas
::
Daxpy
(
m_nquad2
,
input
[
1
+
i
*
totmodes
],
m_base2
.
get
()
+
m_nquad2
,
1
,
&
wsp
[
m_nquad2
*
m_numElmt
]
+
i
*
m_nquad2
,
1
);
// top singular vertex
// (1+c)/2 x (1-b)/2 x (1+a)/2 component
Blas
::
Daxpy
(
m_nquad2
,
input
[
1
+
i
*
totmodes
],
m_base2
.
get
()
+
m_nquad2
,
1
,
&
wsp
[
m_nmodes1
*
m_nquad2
*
m_numElmt
]
+
i
*
m_nquad2
,
1
);
// top singular vertex
// (1+c)/2 x (1+b)/2 x (1+a)/2 component
Blas
::
Daxpy
(
m_nquad2
,
input
[
1
+
i
*
totmodes
],
m_base2
.
get
()
+
m_nquad2
,
1
,
&
wsp
[(
m_nmodes1
+
1
)
*
m_nquad2
*
m_numElmt
]
+
i
*
m_nquad2
,
1
);
}
}
// Perform summation over '1' direction
mode
=
0
;
for
(
i
=
0
;
i
<
m_nmodes0
;
++
i
)
{
Blas
::
Dgemm
(
'N'
,
'T'
,
m_nquad1
,
m_nquad2
*
m_numElmt
,
m_nmodes1
,
1.0
,
m_base1
.
get
(),
m_nquad1
,
wsp
.
get
()
+
mode
*
m_nquad2
*
m_numElmt
,
m_nquad2
*
m_numElmt
,
0.0
,
wsp2
.
get
()
+
i
*
m_nquad1
*
m_nquad2
*
m_numElmt
,
m_nquad1
);
mode
+=
m_nmodes1
;
}
// Perform summation over '0' direction
Blas
::
Dgemm
(
'N'
,
'T'
,
m_nquad0
,
m_nquad1
*
m_nquad2
*
m_numElmt
,
m_nmodes0
,
1.0
,
m_base0
.
get
(),
m_nquad0
,
wsp2
.
get
(),
m_nquad1
*
m_nquad2
*
m_numElmt
,
0.0
,
output
.
get
(),
m_nquad0
);
}
virtual
void
operator
()(
int
dir
,
const
Array
<
OneD
,
const
NekDouble
>
&
input
,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
wsp
)
{
ASSERTL0
(
false
,
"Not valid for this operator."
);
}
protected:
const
int
m_nquad0
;
const
int
m_nquad1
;
const
int
m_nquad2
;
const
int
m_nmodes0
;
const
int
m_nmodes1
;
const
int
m_nmodes2
;
Array
<
OneD
,
const
NekDouble
>
m_base0
;
Array
<
OneD
,
const
NekDouble
>
m_base1
;
Array
<
OneD
,
const
NekDouble
>
m_base2
;
bool
m_sortTopVertex
;
private:
BwdTrans_SumFac_Pyr
(
vector
<
StdRegions
::
StdExpansionSharedPtr
>
pCollExp
,
CoalescedGeomDataSharedPtr
pGeomData
)
:
Operator
(
pCollExp
,
pGeomData
),
m_nquad0
(
m_stdExp
->
GetNumPoints
(
0
)),
m_nquad1
(
m_stdExp
->
GetNumPoints
(
1
)),
m_nquad2
(
m_stdExp
->
GetNumPoints
(
2
)),
m_nmodes0
(
m_stdExp
->
GetBasisNumModes
(
0
)),
m_nmodes1
(
m_stdExp
->
GetBasisNumModes
(
1
)),
m_nmodes2
(
m_stdExp
->
GetBasisNumModes
(
2
)),
m_base0
(
m_stdExp
->
GetBasis
(
0
)
->
GetBdata
()),
m_base1
(
m_stdExp
->
GetBasis
(
1
)
->
GetBdata
()),
m_base2
(
m_stdExp
->
GetBasis
(
2
)
->
GetBdata
())
{
m_wspSize
=
m_numElmt
*
m_nmodes0
*
m_nquad2
*
(
m_nmodes1
+
m_nquad1
);
if
(
m_stdExp
->
GetBasis
(
0
)
->
GetBasisType
()
==
LibUtilities
::
eModified_A
)
{
m_sortTopVertex
=
true
;
}
else
{
m_sortTopVertex
=
false
;
}
}
};
/// Factory initialisation for the BwdTrans_SumFac_Pyr operator
OperatorKey
BwdTrans_SumFac_Pyr
::
m_type
=
GetOperatorFactory
().
RegisterCreatorFunction
(
OperatorKey
(
ePyramid
,
eBwdTrans
,
eSumFac
,
false
),
BwdTrans_SumFac_Pyr
::
create
,
"BwdTrans_SumFac_Pyr"
);
}
}
library/Collections/IProduct.cpp
View file @
4b560009
...
...
@@ -422,6 +422,101 @@ void PrismIProduct(bool sortTopVertex, int numElmt,
}
/**
*
*/
void
PyrIProduct
(
bool
sortTopVertex
,
int
numElmt
,
int
nquad0
,
int
nquad1
,
int
nquad2
,
int
nmodes0
,
int
nmodes1
,
int
nmodes2
,
const
Array
<
OneD
,
const
NekDouble
>
&
base0
,
const
Array
<
OneD
,
const
NekDouble
>
&
base1
,
const
Array
<
OneD
,
const
NekDouble
>
&
base2
,
const
Array
<
OneD
,
const
NekDouble
>
&
jac
,
const
Array
<
OneD
,
const
NekDouble
>
&
input
,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
wsp
)
{
int
totmodes
=
LibUtilities
::
StdPyrData
::
getNumberOfCoefficients
(
nmodes0
,
nmodes1
,
nmodes2
);
int
totpoints
=
nquad0
*
nquad1
*
nquad2
;
int
cnt
;
int
mode
,
mode1
;
ASSERTL1
(
wsp
.
num_elements
()
>=
numElmt
*
(
nquad1
*
nquad2
*
nmodes0
+
nquad2
*
max
(
nquad0
*
nquad1
,
nmodes0
*
nmodes1
)),
"Insufficient workspace size"
);
Vmath
::
Vmul
(
numElmt
*
totpoints
,
jac
,
1
,
input
,
1
,
wsp
,
1
);
Array
<
OneD
,
NekDouble
>
wsp1
=
wsp
+
numElmt
*
nquad2
*
(
max
(
nquad0
*
nquad1
,
nmodes0
*
nmodes1
));
// Perform iproduct with respect to the '0' direction
Blas
::
Dgemm
(
'T'
,
'N'
,
nquad1
*
nquad2
*
numElmt
,
nmodes0
,
nquad0
,
1.0
,
wsp
.
get
(),
nquad0
,
base0
.
get
(),
nquad0
,
0.0
,
wsp1
.
get
(),
nquad1
*
nquad2
*
numElmt
);
// Inner product with respect to the '1' direction
mode
=
0
;
for
(
int
i
=
0
;
i
<
nmodes0
;
++
i
)
{
Blas
::
Dgemm
(
'T'
,
'N'
,
nquad2
*
numElmt
,
nmodes1
,
nquad1
,
1.0
,
wsp1
.
get
()
+
i
*
nquad1
*
nquad2
*
numElmt
,
nquad1
,
base1
.
get
(),
nquad1
,
0.0
,
wsp
.
get
()
+
mode
*
nquad2
*
numElmt
,
nquad2
*
numElmt
);
mode
+=
nmodes1
;
}
// Inner product with respect to the '2' direction
mode
=
mode1
=
cnt
=
0
;
for
(
int
i
=
0
;
i
<
nmodes0
;
++
i
)
{
for
(
int
j
=
0
;
j
<
nmodes1
;
++
j
,
++
cnt
)
{
int
ijmax
=
max
(
i
,
j
);
Blas
::
Dgemm
(
'T'
,
'N'
,
nmodes2
-
ijmax
,
numElmt
,
nquad2
,
1.0
,
base2
.
get
()
+
mode
*
nquad2
,
nquad2
,
wsp
.
get
()
+
cnt
*
nquad2
*
numElmt
,
nquad2
,
0.0
,
output
.
get
()
+
mode1
,
totmodes
);
mode
+=
nmodes2
-
ijmax
;
mode1
+=
nmodes2
-
ijmax
;
}
//increment mode in case order1!=order2
for
(
int
j
=
nmodes1
;
j
<
nmodes2
;
++
j
)
{
int
ijmax
=
max
(
i
,
j
);
mode
+=
nmodes2
-
ijmax
;
}
}
// fix for modified basis for top singular vertex component
// Already have evaluated (1+c)/2 (1-b)/2 (1-a)/2
if
(
sortTopVertex
)
{
for
(
int
n
=
0
;
n
<
numElmt
;
++
n
)
{
// add in (1+c)/2 (1+b)/2 component
output
[
1
+
n
*
totmodes
]
+=
Blas
::
Ddot
(
nquad2
,
base2
.
get
()
+
nquad2
,
1
,
&
wsp
[
nquad2
*
numElmt
+
n
*
nquad2
],
1
);
// add in (1+c)/2 (1-b)/2 (1+a)/2 component
output
[
1
+
n
*
totmodes
]
+=
Blas
::
Ddot
(
nquad2
,
base2
.
get
()
+
nquad2
,
1
,
&
wsp
[
nquad2
*
nmodes1
*
numElmt
+
n
*
nquad2
],
1
);
// add in (1+c)/2 (1+b)/2 (1+a)/2 component
output
[
1
+
n
*
totmodes
]
+=
Blas
::
Ddot
(
nquad2
,
base2
.
get
()
+
nquad2
,
1
,
&
wsp
[
nquad2
*
(
nmodes1
+
1
)
*
numElmt
+
n
*
nquad2
],
1
);
}
}
}
/**
*
*/
...
...
library/Collections/IProduct.h
View file @
4b560009
...
...
@@ -82,6 +82,19 @@ void PrismIProduct(bool sortTopVert, int numElmt,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
wsp
);
void
PyrIProduct
(
bool
sortTopVert
,
int
numElmt
,
int
nquad0
,
int
nquad1
,
int
nquad2
,
int
nmodes0
,
int
nmodes1
,
int
nmodes2
,
const
Array
<
OneD
,
const
NekDouble
>
&
base0
,
const
Array
<
OneD
,
const
NekDouble
>
&
base1
,
const
Array
<
OneD
,
const
NekDouble
>
&
base2
,
const
Array
<
OneD
,
const
NekDouble
>
&
jac
,
const
Array
<
OneD
,
const
NekDouble
>
&
input
,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
wsp
);
void
TetIProduct
(
bool
sortTopEdge
,
int
numElmt
,
int
nquad0
,
int
nquad1
,
int
nquad2
,
int
nmodes0
,
int
nmodes1
,
int
nmodes2
,
...
...
library/Collections/IProductWRTBase.cpp
View file @
4b560009
...
...
@@ -586,7 +586,7 @@ OperatorKey IProductWRTBase_SumFac_Tri::m_type = GetOperatorFactory().
/**
* @brief
Backward transform
operator using sum-factorisation (Hex)
* @brief
Inner Product
operator using sum-factorisation (Hex)
*/
class
IProductWRTBase_SumFac_Hex
:
public
Operator
{
...
...
@@ -765,7 +765,7 @@ OperatorKey IProductWRTBase_SumFac_Tet::m_type = GetOperatorFactory().
/**
* @brief
Backward transform
operator using sum-factorisation (Prism)
* @brief
Inner Product
operator using sum-factorisation (Prism)
*/
class
IProductWRTBase_SumFac_Prism
:
public
Operator
{
...
...
@@ -856,5 +856,99 @@ OperatorKey IProductWRTBase_SumFac_Prism::m_type = GetOperatorFactory().
OperatorKey
(
ePrism
,
eIProductWRTBase
,
eSumFac
,
false
),
IProductWRTBase_SumFac_Prism
::
create
,
"IProductWRTBase_SumFac_Prism"
);
/**
* @brief Inner Product operator using sum-factorisation (Pyr)
*/
class
IProductWRTBase_SumFac_Pyr
:
public
Operator
{
public:
OPERATOR_CREATE
(
IProductWRTBase_SumFac_Pyr
)
virtual
~
IProductWRTBase_SumFac_Pyr
()
{
}
virtual
void
operator
()(
const
Array
<
OneD
,
const
NekDouble
>
&
input
,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
output1
,
Array
<
OneD
,
NekDouble
>
&
output2
,
Array
<
OneD
,
NekDouble
>
&
wsp
)
{
ASSERTL1
(
wsp
.
num_elements
()
==
m_wspSize
,
"Incorrect workspace size"
);
PyrIProduct
(
m_sortTopVertex
,
m_numElmt
,
m_nquad0
,
m_nquad1
,
m_nquad2
,
m_nmodes0
,
m_nmodes1
,
m_nmodes2
,
m_base0
,
m_base1
,
m_base2
,
m_jac
,
input
,
output
,
wsp
);
}
virtual
void
operator
()(
int
dir
,
const
Array
<
OneD
,
const
NekDouble
>
&
input
,
Array
<
OneD
,
NekDouble
>
&
output
,
Array
<
OneD
,
NekDouble
>
&
wsp
)
{
ASSERTL0
(
false
,
"Not valid for this operator."
);
}
protected:
const
int
m_nquad0
;
const
int
m_nquad1
;
const
int
m_nquad2
;
const
int
m_nmodes0
;
const
int
m_nmodes1
;
const
int
m_nmodes2
;
Array
<
OneD
,
const
NekDouble
>
m_jac
;
Array
<
OneD
,
const
NekDouble
>
m_base0
;
Array
<
OneD
,
const
NekDouble
>
m_base1
;
Array
<
OneD
,
const
NekDouble
>
m_base2
;
bool
m_sortTopVertex
;
private:
IProductWRTBase_SumFac_Pyr
(
vector
<
StdRegions
::
StdExpansionSharedPtr
>
pCollExp
,
CoalescedGeomDataSharedPtr
pGeomData
)
:
Operator
(
pCollExp
,
pGeomData
),
m_nquad0
(
m_stdExp
->
GetNumPoints
(
0
)),
m_nquad1
(
m_stdExp
->
GetNumPoints
(
1
)),
m_nquad2
(
m_stdExp
->
GetNumPoints
(
2
)),
m_nmodes0
(
m_stdExp
->
GetBasisNumModes
(
0
)),
m_nmodes1
(
m_stdExp
->
GetBasisNumModes
(
1
)),
m_nmodes2
(
m_stdExp
->
GetBasisNumModes
(
2
)),
m_base0
(
m_stdExp
->
GetBasis
(
0
)
->
GetBdata
()),
m_base1
(
m_stdExp
->
GetBasis
(
1
)
->
GetBdata
()),
m_base2
(
m_stdExp
->
GetBasis
(
2
)
->
GetBdata
())
{
m_jac
=
pGeomData
->
GetJacWithStdWeights
(
pCollExp
);
m_wspSize
=
m_numElmt
*
m_nquad2
*
(
max
(
m_nquad0
*
m_nquad1
,
m_nmodes0
*
m_nmodes1
))
+
m_nquad1
*
m_nquad2
*
m_numElmt
*
m_nmodes0
;
if
(
m_stdExp
->
GetBasis
(
0
)
->
GetBasisType
()
==
LibUtilities
::
eModified_A
)
{
m_sortTopVertex
=
true
;
}
else
{
m_sortTopVertex
=
false
;
}
}
};
/// Factory initialisation for the IProductWRTBase_SumFac_Pyr operator
OperatorKey
IProductWRTBase_SumFac_Pyr
::
m_type
=
GetOperatorFactory
().
RegisterCreatorFunction
(
OperatorKey
(
ePyramid
,
eIProductWRTBase
,
eSumFac
,
false
),
IProductWRTBase_SumFac_Pyr
::
create
,
"IProductWRTBase_SumFac_Pyr"
);
}
}
library/Collections/IProductWRTDerivBase.cpp
View file @
4b560009
This diff is collapsed.
Click to expand it.
library/Collections/PhysDeriv.cpp
View file @
4b560009
...
...
@@ -938,7 +938,7 @@ class PhysDeriv_SumFac_Hex : public Operator
Array
<
OneD
,
Array
<
OneD
,
NekDouble
>
>
out
(
3
);
out
[
0
]
=
output0
;
out
[
1
]
=
output1
;
out
[
2
]
=
output2
;
for
(
int
i
=
0
;
i
<
m_dim
;
++
i
)
for
(
int
i
=
0
;
i
<
3
;
++
i
)
{
Diff
[
i
]
=
wsp
+
i
*
ntot
;
}
...
...
@@ -967,10 +967,10 @@ class PhysDeriv_SumFac_Hex : public Operator
// calculate full derivative
for
(
int
i
=
0
;
i
<
m_coordim
;
++
i
)
{
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
i
*
m_dim
],
1
,
Diff
[
0
],
1
,
out
[
i
],
1
);
for
(
int
j
=
1
;
j
<
m_dim
;
++
j
)
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
i
*
3
],
1
,
Diff
[
0
],
1
,
out
[
i
],
1
);
for
(
int
j
=
1
;
j
<
3
;
++
j
)
{
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
i
*
m_dim
+
j
],
1
,
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
i
*
3
+
j
],
1
,
Diff
[
j
],
1
,
out
[
i
],
1
,
out
[
i
],
1
);
...
...
@@ -989,7 +989,7 @@ class PhysDeriv_SumFac_Hex : public Operator
Array
<
OneD
,
NekDouble
>
tmp0
,
tmp1
,
tmp2
;
Array
<
OneD
,
Array
<
OneD
,
NekDouble
>
>
Diff
(
3
);
for
(
int
i
=
0
;
i
<
m_dim
;
++
i
)
for
(
int
i
=
0
;
i
<
3
;
++
i
)
{
Diff
[
i
]
=
wsp
+
i
*
ntot
;
}
...
...
@@ -1016,10 +1016,10 @@ class PhysDeriv_SumFac_Hex : public Operator
}
// calculate full derivative
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
dir
*
m_dim
],
1
,
Diff
[
0
],
1
,
output
,
1
);
for
(
int
j
=
1
;
j
<
m_dim
;
++
j
)
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
dir
*
3
],
1
,
Diff
[
0
],
1
,
output
,
1
);
for
(
int
j
=
1
;
j
<
3
;
++
j
)
{
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
dir
*
m_dim
+
j
],
1
,
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
dir
*
3
+
j
],
1
,
Diff
[
j
],
1
,
output
,
1
,
output
,
1
);
...
...
@@ -1028,7 +1028,6 @@ class PhysDeriv_SumFac_Hex : public Operator
protected:
Array
<
TwoD
,
const
NekDouble
>
m_derivFac
;
int
m_dim
;
int
m_coordim
;
const
int
m_nquad0
;
const
int
m_nquad1
;
...
...
@@ -1048,7 +1047,6 @@ class PhysDeriv_SumFac_Hex : public Operator
{
LibUtilities
::
PointsKeyVector
PtsKey
=
m_stdExp
->
GetPointsKeys
();
m_dim
=
PtsKey
.
size
();
m_coordim
=
m_stdExp
->
GetCoordim
();
m_derivFac
=
pGeomData
->
GetDerivFactors
(
pCollExp
);
...
...
@@ -1096,7 +1094,7 @@ class PhysDeriv_SumFac_Tet : public Operator
Array
<
OneD
,
Array
<
OneD
,
NekDouble
>
>
out
(
3
);
out
[
0
]
=
output0
;
out
[
1
]
=
output1
;
out
[
2
]
=
output2
;
for
(
int
i
=
0
;
i
<
m_dim
;
++
i
)
for
(
int
i
=
0
;
i
<
3
;
++
i
)
{
Diff
[
i
]
=
wsp
+
i
*
ntot
;
}
...
...
@@ -1161,10 +1159,10 @@ class PhysDeriv_SumFac_Tet : public Operator
// calculate full derivative
for
(
int
i
=
0
;
i
<
m_coordim
;
++
i
)
{
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
i
*
m_dim
],
1
,
Diff
[
0
],
1
,
out
[
i
],
1
);
for
(
int
j
=
1
;
j
<
m_dim
;
++
j
)
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
i
*
3
],
1
,
Diff
[
0
],
1
,
out
[
i
],
1
);
for
(
int
j
=
1
;
j
<
3
;
++
j
)
{
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
i
*
m_dim
+
j
],
1
,
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
i
*
3
+
j
],
1
,
Diff
[
j
],
1
,
out
[
i
],
1
,
out
[
i
],
1
);
}
}
...
...
@@ -1181,7 +1179,7 @@ class PhysDeriv_SumFac_Tet : public Operator
Array
<
OneD
,
NekDouble
>
tmp0
,
tmp1
,
tmp2
;
Array
<
OneD
,
Array
<
OneD
,
NekDouble
>
>
Diff
(
3
);
for
(
int
i
=
0
;
i
<
m_dim
;
++
i
)
for
(
int
i
=
0
;
i
<
3
;
++
i
)
{
Diff
[
i
]
=
wsp
+
i
*
ntot
;
}
...
...
@@ -1244,17 +1242,16 @@ class PhysDeriv_SumFac_Tet : public Operator
}
// calculate full derivative
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
dir
*
m_dim
],
1
,
Diff
[
0
],
1
,
output
,
1
);
for
(
int
j
=
1
;
j
<
m_dim
;
++
j
)
Vmath
::
Vmul
(
ntot
,
m_derivFac
[
dir
*
3
],
1
,
Diff
[
0
],
1
,
output
,
1
);
for
(
int
j
=
1
;
j
<
3
;
++
j
)
{
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
dir
*
m_dim
+
j
],
1
,
Vmath
::
Vvtvp
(
ntot
,
m_derivFac
[
dir
*
3
+
j
],
1
,
Diff
[
j
],
1
,
output
,
1
,
output
,
1
);
}
}
protected:
Array
<
TwoD
,
const
NekDouble
>
m_derivFac
;
int
m_dim
;
int
m_coordim
;
const
int
m_nquad0
;
const
int
m_nquad1
;
...
...
@@ -1278,7 +1275,6 @@ class PhysDeriv_SumFac_Tet : public Operator
{
LibUtilities
::
PointsKeyVector
PtsKey
=
m_stdExp
->
GetPointsKeys
();
m_dim
=
PtsKey
.
size
();
m_coordim
=
m_stdExp
->
GetCoordim
();
m_derivFac
=
pGeomData
->
GetDerivFactors
(
pCollExp
);
...
...
@@ -1358,7 +1354,7 @@ class PhysDeriv_SumFac_Prism : public Operator
Array
<
OneD
,
Array
<
OneD
,
NekDouble
>
>
out
(
3
);
out
[
0
]
=
output0
;
out
[
1
]
=
output1
;
out
[
2
]
=
output2
;
for
(
int
i
=
0
;
i
<
m_dim
;
++
i
)
for
(
int
i
=
0
;
i
<
3
;
++
i
)