Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in / Register
Toggle navigation
Nektar
Project overview
Project overview
Details
Activity
Releases
Repository
Repository
Files
Commits
Branches
Tags
Contributors
Graph
Compare
Locked Files
Issues
0
Issues
0
List
Boards
Labels
Service Desk
Milestones
Iterations
Merge Requests
0
Merge Requests
0
Requirements
Requirements
List
Security & Compliance
Security & Compliance
Dependency List
License Compliance
Operations
Operations
Incidents
Packages & Registries
Packages & Registries
Container Registry
Analytics
Analytics
Code Review
Insights
Issue
Repository
Value Stream
Wiki
Wiki
Members
Members
Collapse sidebar
Close sidebar
Activity
Graph
Create a new issue
Commits
Issue Boards
Open sidebar
Jennifer Ryan
Nektar
Commits
4c83a458
Commit
4c83a458
authored
Apr 10, 2017
by
Dave Moxey
Browse files
Options
Browse Files
Download
Plain Diff
Merge branch 'feature/SumFacPyr' into 'master'
Feature/sum fac pyr See merge request !750
parents
236cb5f9
9d6a3423
Changes
73
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
73 changed files
with
4557 additions
and
1032 deletions
+4557
-1032
CHANGELOG.md
CHANGELOG.md
+4
-0
docs/doxygen/Doxyfile.in
docs/doxygen/Doxyfile.in
+1
-0
library/Collections/BwdTrans.cpp
library/Collections/BwdTrans.cpp
+164
-1
library/Collections/IProduct.cpp
library/Collections/IProduct.cpp
+95
-0
library/Collections/IProduct.h
library/Collections/IProduct.h
+13
-0
library/Collections/IProductWRTBase.cpp
library/Collections/IProductWRTBase.cpp
+96
-2
library/Collections/IProductWRTDerivBase.cpp
library/Collections/IProductWRTDerivBase.cpp
+383
-85
library/Collections/PhysDeriv.cpp
library/Collections/PhysDeriv.cpp
+247
-30
library/Demos/LocalRegions/LocProject1D.cpp
library/Demos/LocalRegions/LocProject1D.cpp
+8
-8
library/Demos/LocalRegions/LocProject2D.cpp
library/Demos/LocalRegions/LocProject2D.cpp
+13
-12
library/Demos/LocalRegions/LocProject3D.cpp
library/Demos/LocalRegions/LocProject3D.cpp
+18
-13
library/Demos/LocalRegions/LocProject_Diff1D.cpp
library/Demos/LocalRegions/LocProject_Diff1D.cpp
+9
-9
library/Demos/LocalRegions/LocProject_Diff2D.cpp
library/Demos/LocalRegions/LocProject_Diff2D.cpp
+13
-14
library/Demos/LocalRegions/LocProject_Diff3D.cpp
library/Demos/LocalRegions/LocProject_Diff3D.cpp
+18
-10
library/Demos/LocalRegions/Tests/LocProject2D_Def_Quad_Lagrange_Basis_P6_Q7.tst
...ions/Tests/LocProject2D_Def_Quad_Lagrange_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject2D_Quad_Lagrange_Basis_P6_Q7.tst
...lRegions/Tests/LocProject2D_Quad_Lagrange_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject2D_Tri_Nodal_Basis_P6_Q7.tst
...LocalRegions/Tests/LocProject2D_Tri_Nodal_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject3D_Def_Hex_Lagrange_Basis_P6_Q7.tst
...gions/Tests/LocProject3D_Def_Hex_Lagrange_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject3D_Hex_Lagrange_Basis_P6_Q7.tst
...alRegions/Tests/LocProject3D_Hex_Lagrange_Basis_P6_Q7.tst
+3
-3
library/Demos/LocalRegions/Tests/LocProject3D_Pyr_Mod_Basis_P6_Q7.tst
...s/LocalRegions/Tests/LocProject3D_Pyr_Mod_Basis_P6_Q7.tst
+3
-3
library/Demos/LocalRegions/Tests/LocProject_Diff2D_Lin_Deformed_Quad_Lagrange_Basis_P6_Q7.tst
...Project_Diff2D_Lin_Deformed_Quad_Lagrange_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject_Diff2D_Reg_Quad_Lagrange_Basis_P6_Q=7.tst
...ests/LocProject_Diff2D_Reg_Quad_Lagrange_Basis_P6_Q=7.tst
+2
-2
library/Demos/LocalRegions/Tests/LocProject_Diff2D_Tri_Nodal_Basis_P6_Q7.tst
...Regions/Tests/LocProject_Diff2D_Tri_Nodal_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject_Diff3D_Lin_Deformed_Hex_Lagrange_Basis_P6_Q7.tst
...cProject_Diff3D_Lin_Deformed_Hex_Lagrange_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject_Diff3D_Reg_Hex_Lagrange_Basis_P6_Q7.tst
.../Tests/LocProject_Diff3D_Reg_Hex_Lagrange_Basis_P6_Q7.tst
+1
-1
library/Demos/LocalRegions/Tests/LocProject_Diff3D_Reg_Pyr_Mod_Basis_P6_Q7.tst
...gions/Tests/LocProject_Diff3D_Reg_Pyr_Mod_Basis_P6_Q7.tst
+3
-3
library/Demos/MultiRegions/Tests/Helmholtz3D_CG_Pyr.tst
library/Demos/MultiRegions/Tests/Helmholtz3D_CG_Pyr.tst
+2
-2
library/Demos/MultiRegions/Tests/Helmholtz3D_CG_Pyr_Deformed.tst
.../Demos/MultiRegions/Tests/Helmholtz3D_CG_Pyr_Deformed.tst
+2
-2
library/Demos/StdRegions/CMakeLists.txt
library/Demos/StdRegions/CMakeLists.txt
+1
-0
library/Demos/StdRegions/StdProject1D.cpp
library/Demos/StdRegions/StdProject1D.cpp
+8
-8
library/Demos/StdRegions/StdProject2D.cpp
library/Demos/StdRegions/StdProject2D.cpp
+13
-12
library/Demos/StdRegions/StdProject3D.cpp
library/Demos/StdRegions/StdProject3D.cpp
+55
-32
library/Demos/StdRegions/StdProject_Diff1D.cpp
library/Demos/StdRegions/StdProject_Diff1D.cpp
+8
-8
library/Demos/StdRegions/StdProject_Diff2D.cpp
library/Demos/StdRegions/StdProject_Diff2D.cpp
+11
-11
library/Demos/StdRegions/StdProject_Diff3D.cpp
library/Demos/StdRegions/StdProject_Diff3D.cpp
+25
-12
library/Demos/StdRegions/Tests/StdProject2D_Quad_Fourier_P6_Q8.tst
...emos/StdRegions/Tests/StdProject2D_Quad_Fourier_P6_Q8.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject2D_Quad_Fourier_Single_Mode_P2_Q6.tst
...ons/Tests/StdProject2D_Quad_Fourier_Single_Mode_P2_Q6.tst
+2
-2
library/Demos/StdRegions/Tests/StdProject2D_Quad_Lagrange_P6_Q7.tst
...mos/StdRegions/Tests/StdProject2D_Quad_Lagrange_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject2D_Tri_Nodal_P6_Q7.tst
...y/Demos/StdRegions/Tests/StdProject2D_Tri_Nodal_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject3D_Hex_Chebyshev_P6_Q7.tst
...mos/StdRegions/Tests/StdProject3D_Hex_Chebyshev_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject3D_Hex_Lagrange_P6_Q7.tst
...emos/StdRegions/Tests/StdProject3D_Hex_Lagrange_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject3D_Hex_Legendre_P6_Q7.tst
...emos/StdRegions/Tests/StdProject3D_Hex_Legendre_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject3D_Pyr_Mod_P6_Q7.tst
...ary/Demos/StdRegions/Tests/StdProject3D_Pyr_Mod_P6_Q7.tst
+3
-3
library/Demos/StdRegions/Tests/StdProject3D_Pyr_Ortho_P6_Q7.tst
...y/Demos/StdRegions/Tests/StdProject3D_Pyr_Ortho_P6_Q7.tst
+14
-0
library/Demos/StdRegions/Tests/StdProject_Diff2D_Quad_Fourier_P6_Q8.tst
...StdRegions/Tests/StdProject_Diff2D_Quad_Fourier_P6_Q8.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject_Diff2D_Quad_Lagrange_P6_Q7.tst
...tdRegions/Tests/StdProject_Diff2D_Quad_Lagrange_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject_Diff2D_Tri_Nodal_P6_Q7.tst
...os/StdRegions/Tests/StdProject_Diff2D_Tri_Nodal_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject_Diff3D_Hex_Chebyshev_P6_Q7.tst
...tdRegions/Tests/StdProject_Diff3D_Hex_Chebyshev_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject_Diff3D_Hex_Lagrange_P6_Q7.tst
...StdRegions/Tests/StdProject_Diff3D_Hex_Lagrange_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject_Diff3D_Hex_Legendre_P6_Q7.tst
...StdRegions/Tests/StdProject_Diff3D_Hex_Legendre_P6_Q7.tst
+1
-1
library/Demos/StdRegions/Tests/StdProject_Diff3D_Pyr_Mod_P6_Q7.tst
...emos/StdRegions/Tests/StdProject_Diff3D_Pyr_Mod_P6_Q7.tst
+1
-1
library/LibUtilities/BasicUtils/ShapeType.hpp
library/LibUtilities/BasicUtils/ShapeType.hpp
+6
-11
library/LibUtilities/Foundations/Basis.cpp
library/LibUtilities/Foundations/Basis.cpp
+175
-9
library/LibUtilities/Foundations/Basis.h
library/LibUtilities/Foundations/Basis.h
+4
-1
library/LibUtilities/Foundations/BasisType.h
library/LibUtilities/Foundations/BasisType.h
+2
-0
library/LibUtilities/Foundations/Foundations.hpp
library/LibUtilities/Foundations/Foundations.hpp
+2
-0
library/LocalRegions/PyrExp.cpp
library/LocalRegions/PyrExp.cpp
+217
-13
library/LocalRegions/PyrExp.h
library/LocalRegions/PyrExp.h
+18
-0
library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
+3
-2
library/SpatialDomains/MeshGraph.cpp
library/SpatialDomains/MeshGraph.cpp
+2
-2
library/SpatialDomains/PyrGeom.cpp
library/SpatialDomains/PyrGeom.cpp
+3
-3
library/StdRegions/StdExpansion.cpp
library/StdRegions/StdExpansion.cpp
+1
-1
library/StdRegions/StdExpansion3D.cpp
library/StdRegions/StdExpansion3D.cpp
+2
-0
library/StdRegions/StdPrismExp.cpp
library/StdRegions/StdPrismExp.cpp
+4
-4
library/StdRegions/StdPyrExp.cpp
library/StdRegions/StdPyrExp.cpp
+676
-668
library/StdRegions/StdPyrExp.h
library/StdRegions/StdPyrExp.h
+11
-4
library/UnitTests/Collections/CMakeLists.txt
library/UnitTests/Collections/CMakeLists.txt
+1
-0
library/UnitTests/Collections/TestPyrCollection.cpp
library/UnitTests/Collections/TestPyrCollection.cpp
+1973
-0
solvers/IncNavierStokesSolver/CMakeLists.txt
solvers/IncNavierStokesSolver/CMakeLists.txt
+1
-0
solvers/IncNavierStokesSolver/Tests/ChanFlow_m8_BodyForce.xml
...ers/IncNavierStokesSolver/Tests/ChanFlow_m8_BodyForce.xml
+1
-0
solvers/IncNavierStokesSolver/Tests/ChanFlow_m8_BodyForce_par.tst
...IncNavierStokesSolver/Tests/ChanFlow_m8_BodyForce_par.tst
+6
-6
solvers/IncNavierStokesSolver/Tests/Pyr_channel_SVV.tst
solvers/IncNavierStokesSolver/Tests/Pyr_channel_SVV.tst
+24
-0
solvers/IncNavierStokesSolver/Tests/Pyr_channel_SVV.xml
solvers/IncNavierStokesSolver/Tests/Pyr_channel_SVV.xml
+167
-0
No files found.
CHANGELOG.md
View file @
4c83a458
...
...
@@ -7,6 +7,10 @@ v4.5.0
-
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)
...
...
docs/doxygen/Doxyfile.in
View file @
4c83a458
...
...
@@ -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 @
4c83a458
...
...
@@ -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 @
4c83a458
...
...
@@ -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 @
4c83a458
...
...
@@ -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 @
4c83a458
...
...
@@ -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 @
4c83a458
This diff is collapsed.
Click to expand it.
library/Collections/PhysDeriv.cpp
View file @
4c83a458
...
...
@@ -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
)