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
612a6e36
Commit
612a6e36
authored
Mar 04, 2013
by
Chris Cantwell
Browse files
Merge branch 'feature/FieldIO' of localhost:nektar
parents
d3948e5b
26575e56
Changes
235
Hide whitespace changes
Inline
Side-by-side
library/Demos/LocalRegions/LocPrismExpDemo.cpp
View file @
612a6e36
...
...
@@ -55,7 +55,7 @@ int main(int argc, char *argv[])
exit
(
1
);
}
StdRegions
::
Expansion
Type
regionShape
=
StdRegion
s
::
ePrism
;
LibUtilities
::
Shape
Type
regionShape
=
LibUtilitie
s
::
ePrism
;
int
bType_x_val
=
atoi
(
argv
[
1
]);
int
bType_y_val
=
atoi
(
argv
[
2
]);
int
bType_z_val
=
atoi
(
argv
[
3
]);
...
...
@@ -72,7 +72,7 @@ int main(int argc, char *argv[])
}
// Check to see that correct Expansions are used
if
(
regionShape
==
StdRegion
s
::
ePrism
)
if
(
regionShape
==
LibUtilitie
s
::
ePrism
)
{
if
(
(
bType_x
==
LibUtilities
::
eOrtho_B
)
||
(
bType_x
==
LibUtilities
::
eModified_B
)
)
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 1 cannot be of type Ortho_B or Modified_B"
);
...
...
@@ -105,7 +105,7 @@ int main(int argc, char *argv[])
// Define a 3D expansion based on basis definition
StdRegions
::
StdExpansion3D
*
lpr
=
0
;
if
(
regionShape
==
StdRegion
s
::
ePrism
)
if
(
regionShape
==
LibUtilitie
s
::
ePrism
)
{
// //////////////////////////////////////////////////////
// Set up Prism vertex coordinates
...
...
@@ -197,7 +197,7 @@ int main(int argc, char *argv[])
Array
<
OneD
,
StdRegions
::
StdExpansion3DSharedPtr
>
xMap
(
3
);
for
(
int
i
=
0
;
i
<
3
;
++
i
){
xMap
[
i
]
=
MemoryManager
<
StdRegions
::
StdPrismExp
>::
AllocateSharedPtr
(
basisKey_x
,
basisKey_y
,
basisKey_z
);
xMap
[
i
]
=
MemoryManager
<
StdRegions
::
StdPrismExp
>::
AllocateSharedPtr
(
basisKey_x
,
basisKey_y
,
basisKey_z
);
}
...
...
@@ -250,7 +250,8 @@ int main(int argc, char *argv[])
t
[
1
]
=
0.5
;
t
[
2
]
=
0.2
;
if
(
regionShape
==
StdRegions
::
ePrism
)
{
if
(
regionShape
==
LibUtilities
::
ePrism
)
{
solution
[
0
]
=
Prism_sol
(
t
[
0
],
t
[
1
],
t
[
2
],
P
,
Q
,
R
,
bType_x
,
bType_y
,
bType_z
);
}
...
...
library/Demos/LocalRegions/LocProject2D.cpp
View file @
612a6e36
...
...
@@ -25,7 +25,7 @@ int main(int argc, char *argv[])
LibUtilities
::
PointsType
Qtype1
,
Qtype2
;
LibUtilities
::
BasisType
btype1
,
btype2
;
LibUtilities
::
PointsType
NodalType
;
StdRegions
::
Expansion
Type
regionshape
;
LibUtilities
::
Shape
Type
regionshape
;
StdRegions
::
StdExpansion2D
*
E
;
Array
<
OneD
,
NekDouble
>
sol
;
Array
<
OneD
,
NekDouble
>
coords
(
8
);
...
...
@@ -41,8 +41,8 @@ int main(int argc, char *argv[])
fprintf
(
stderr
,
"Where RegionShape is an integer value which "
"dictates the region shape:
\n
"
);
fprintf
(
stderr
,
"
\t
Triangle =
2
\n
"
);
fprintf
(
stderr
,
"
\t
Quadrilateral =
3
\n
"
);
fprintf
(
stderr
,
"
\t
Triangle =
3
\n
"
);
fprintf
(
stderr
,
"
\t
Quadrilateral =
4
\n
"
);
fprintf
(
stderr
,
"Where type is an integer value which "
"dictates the basis as:
\n
"
);
...
...
@@ -64,10 +64,10 @@ int main(int argc, char *argv[])
exit
(
1
);
}
regionshape
=
(
StdRegions
::
Expansion
Type
)
atoi
(
argv
[
1
]);
regionshape
=
(
LibUtilities
::
Shape
Type
)
atoi
(
argv
[
1
]);
// Check to see if 2D region
if
((
regionshape
!=
StdRegion
s
::
eTriangle
)
&&
(
regionshape
!=
StdRegion
s
::
eQuadrilateral
))
if
((
regionshape
!=
LibUtilitie
s
::
eTriangle
)
&&
(
regionshape
!=
LibUtilitie
s
::
eQuadrilateral
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"This shape is not a 2D region"
);
}
...
...
@@ -99,7 +99,7 @@ int main(int argc, char *argv[])
// Check to see that correct Expansions are used
switch
(
regionshape
)
{
case
StdRegion
s
::
eTriangle
:
case
LibUtilitie
s
::
eTriangle
:
if
((
btype1
==
LibUtilities
::
eOrtho_B
)
||
(
btype1
==
LibUtilities
::
eModified_B
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
...
...
@@ -107,7 +107,7 @@ int main(int argc, char *argv[])
}
break
;
case
StdRegion
s
::
eQuadrilateral
:
case
LibUtilitie
s
::
eQuadrilateral
:
if
((
btype1
==
LibUtilities
::
eOrtho_B
)
||
(
btype1
==
LibUtilities
::
eOrtho_C
)
||
(
btype1
==
LibUtilities
::
eModified_B
)
||
(
btype1
==
LibUtilities
::
eModified_C
))
{
...
...
@@ -145,7 +145,7 @@ int main(int argc, char *argv[])
if
(
btype2
!=
LibUtilities
::
eFourier
)
{
if
(
regionshape
==
StdRegion
s
::
eTriangle
)
{
if
(
regionshape
==
LibUtilitie
s
::
eTriangle
)
{
Qtype2
=
LibUtilities
::
eGaussRadauMAlpha1Beta0
;
}
else
...
...
@@ -163,9 +163,8 @@ int main(int argc, char *argv[])
switch
(
regionshape
)
{
case
StdRegion
s
::
eTriangle
:
case
LibUtilitie
s
::
eTriangle
:
{
coords
[
0
]
=
atof
(
argv
[
8
]);
coords
[
1
]
=
atof
(
argv
[
9
]);
coords
[
2
]
=
atof
(
argv
[
10
]);
...
...
@@ -225,7 +224,7 @@ int main(int argc, char *argv[])
}
break
;
case
StdRegion
s
::
eQuadrilateral
:
case
LibUtilitie
s
::
eQuadrilateral
:
{
// Gather coordinates
coords
[
0
]
=
atof
(
argv
[
8
]);
...
...
@@ -319,7 +318,7 @@ int main(int argc, char *argv[])
x
[
0
]
=
(
coords
[
0
]
+
coords
[
2
])
*
0.5
;
x
[
1
]
=
(
coords
[
1
]
+
coords
[
5
])
*
0.5
;
if
(
regionshape
==
StdRegion
s
::
eTriangle
)
if
(
regionshape
==
LibUtilitie
s
::
eTriangle
)
{
sol
[
0
]
=
Tri_sol
(
x
[
0
],
x
[
1
],
order1
,
order2
);
}
...
...
library/Demos/LocalRegions/LocProject3D.cpp
View file @
612a6e36
...
...
@@ -58,7 +58,7 @@ int main(int argc, char *argv[]){
PointsType
Qtype1
,
Qtype2
,
Qtype3
;
BasisType
btype1
,
btype2
,
btype3
;
Expansion
Type
regionshape
;
Shape
Type
regionshape
;
StdExpansion
*
E
;
Array
<
OneD
,
NekDouble
>
sol
;
...
...
@@ -69,9 +69,9 @@ int main(int argc, char *argv[]){
"x3 y3 z3 [x4 y4 z4...]
\n
"
);
fprintf
(
stderr
,
"Where RegionShape is an integer value which "
"dictates the region shape:
\n
"
);
fprintf
(
stderr
,
"
\t
Tetrahedron =
4
\n
"
);
fprintf
(
stderr
,
"
\t
Prism =
6
\n
"
);
fprintf
(
stderr
,
"
\t
Hexahedron =
7
\n
"
);
fprintf
(
stderr
,
"
\t
Tetrahedron =
5
\n
"
);
fprintf
(
stderr
,
"
\t
Prism =
7
\n
"
);
fprintf
(
stderr
,
"
\t
Hexahedron =
8
\n
"
);
fprintf
(
stderr
,
"Where type is an integer value which "
...
...
@@ -91,12 +91,12 @@ int main(int argc, char *argv[]){
exit
(
1
);
}
regionshape
=
(
StdRegions
::
Expansion
Type
)
atoi
(
argv
[
1
]);
regionshape
=
(
LibUtilities
::
Shape
Type
)
atoi
(
argv
[
1
]);
// Check to see if 3D region
if
(
regionshape
!=
StdRegion
s
::
eTetrahedron
&&
regionshape
!=
StdRegion
s
::
ePrism
&&
regionshape
!=
StdRegion
s
::
eHexahedron
)
if
(
regionshape
!=
LibUtilitie
s
::
eTetrahedron
&&
regionshape
!=
LibUtilitie
s
::
ePrism
&&
regionshape
!=
LibUtilitie
s
::
eHexahedron
)
{
NEKERROR
(
ErrorUtil
::
efatal
,
"This shape is not a 3D region"
);
}
...
...
@@ -111,7 +111,7 @@ int main(int argc, char *argv[]){
// Check to see that correct Expansions are used
switch
(
regionshape
)
{
case
StdRegion
s
::
eTetrahedron
:
case
LibUtilitie
s
::
eTetrahedron
:
if
((
btype1
==
eOrtho_B
)
||
(
btype1
==
eOrtho_C
)
||
(
btype1
==
eModified_B
)
||
(
btype1
==
eModified_C
))
{
...
...
@@ -134,7 +134,7 @@ int main(int argc, char *argv[]){
"or Modified_B"
);
}
break
;
case
StdRegion
s
::
ePrism
:
case
LibUtilitie
s
::
ePrism
:
if
((
btype1
==
eOrtho_B
)
||
(
btype1
==
eOrtho_C
)
||
(
btype1
==
eModified_B
)
||
(
btype1
==
eModified_C
))
{
...
...
@@ -157,7 +157,7 @@ int main(int argc, char *argv[]){
"or Modified_C"
);
}
break
;
case
StdRegion
s
::
eHexahedron
:
case
LibUtilitie
s
::
eHexahedron
:
if
((
btype1
==
eOrtho_B
)
||
(
btype1
==
eOrtho_C
)
||
(
btype1
==
eModified_B
)
||
(
btype1
==
eModified_C
))
{
...
...
@@ -204,7 +204,7 @@ int main(int argc, char *argv[]){
if
(
btype2
!=
LibUtilities
::
eFourier
)
{
if
(
regionshape
==
StdRegion
s
::
eTetrahedron
)
{
if
(
regionshape
==
LibUtilitie
s
::
eTetrahedron
)
{
Qtype2
=
LibUtilities
::
eGaussRadauMAlpha1Beta0
;
}
else
...
...
@@ -219,10 +219,11 @@ int main(int argc, char *argv[]){
if
(
btype3
!=
LibUtilities
::
eFourier
)
{
if
(
regionshape
==
StdRegions
::
eTetrahedron
)
{
if
(
regionshape
==
LibUtilities
::
eTetrahedron
)
{
Qtype3
=
LibUtilities
::
eGaussRadauMAlpha2Beta0
;
}
else
if
(
regionshape
==
StdRegion
s
::
ePrism
)
else
if
(
regionshape
==
LibUtilitie
s
::
ePrism
)
{
Qtype3
=
LibUtilities
::
eGaussRadauMAlpha1Beta0
;
}
...
...
@@ -251,7 +252,7 @@ int main(int argc, char *argv[]){
switch
(
regionshape
)
{
case
StdRegion
s
::
eTetrahedron
:
case
LibUtilitie
s
::
eTetrahedron
:
{
SpatialDomains
::
TetGeomSharedPtr
geom
=
CreateTetGeom
(
argc
,
argv
);
E
=
new
LocalRegions
::
TetExp
(
Bkey1
,
Bkey2
,
Bkey3
,
geom
);
...
...
@@ -266,7 +267,7 @@ int main(int argc, char *argv[]){
//----------------------------------------------
}
break
;
case
StdRegion
s
::
ePrism
:
case
LibUtilitie
s
::
ePrism
:
{
SpatialDomains
::
PrismGeomSharedPtr
geom
=
CreatePrismGeom
(
argc
,
argv
);
E
=
new
LocalRegions
::
PrismExp
(
Bkey1
,
Bkey2
,
Bkey3
,
geom
);
...
...
@@ -281,7 +282,7 @@ int main(int argc, char *argv[]){
//----------------------------------------------
}
break
;
case
StdRegion
s
::
eHexahedron
:
case
LibUtilitie
s
::
eHexahedron
:
{
SpatialDomains
::
HexGeomSharedPtr
geom
=
CreateHexGeom
(
argc
,
argv
);
E
=
new
LocalRegions
::
HexExp
(
Bkey1
,
Bkey2
,
Bkey3
,
geom
);
...
...
library/Demos/LocalRegions/LocProject_Diff2D.cpp
View file @
612a6e36
...
...
@@ -35,7 +35,7 @@ int main(int argc, char *argv[])
LibUtilities
::
PointsType
Qtype1
,
Qtype2
;
LibUtilities
::
BasisType
btype1
,
btype2
;
LibUtilities
::
PointsType
NodalType
;
StdRegions
::
Expansion
Type
regionshape
;
LibUtilities
::
Shape
Type
regionshape
;
StdRegions
::
StdExpansion2D
*
E
;
Array
<
OneD
,
NekDouble
>
sol
,
x
,
y
,
dx
,
dy
;
Array
<
OneD
,
NekDouble
>
coords
(
8
);
...
...
@@ -51,8 +51,8 @@ int main(int argc, char *argv[])
fprintf
(
stderr
,
"Where RegionShape is an integer value which "
"dictates the region shape:
\n
"
);
fprintf
(
stderr
,
"
\t
Triangle =
2
\n
"
);
fprintf
(
stderr
,
"
\t
Quadrilateral =
3
\n
"
);
fprintf
(
stderr
,
"
\t
Triangle =
3
\n
"
);
fprintf
(
stderr
,
"
\t
Quadrilateral =
4
\n
"
);
fprintf
(
stderr
,
"Where type is an integer value which "
"dictates the basis as:
\n
"
);
...
...
@@ -74,10 +74,11 @@ int main(int argc, char *argv[])
exit
(
1
);
}
regionshape
=
(
StdRegions
::
Expansion
Type
)
atoi
(
argv
[
1
]);
regionshape
=
(
LibUtilities
::
Shape
Type
)
atoi
(
argv
[
1
]);
// Check to see if 2D region
if
((
regionshape
!=
StdRegions
::
eTriangle
)
&&
(
regionshape
!=
StdRegions
::
eQuadrilateral
))
if
((
regionshape
!=
LibUtilities
::
eTriangle
)
&&
(
regionshape
!=
LibUtilities
::
eQuadrilateral
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"This shape is not a 2D region"
);
}
...
...
@@ -110,35 +111,35 @@ int main(int argc, char *argv[])
// Check to see that correct Expansions are used
switch
(
regionshape
)
{
case
StdRegion
s
::
eTriangle
:
if
((
btype1
==
LibUtilities
::
eOrtho_B
)
||
(
btype1
==
LibUtilities
::
eModified_B
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 1 cannot be of type Ortho_B or Modified_B"
);
}
case
LibUtilitie
s
::
eTriangle
:
if
((
btype1
==
LibUtilities
::
eOrtho_B
)
||
(
btype1
==
LibUtilities
::
eModified_B
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 1 cannot be of type Ortho_B or Modified_B"
);
}
break
;
case
StdRegion
s
::
eQuadrilateral
:
if
((
btype1
==
LibUtilities
::
eOrtho_B
)
||
(
btype1
==
LibUtilities
::
eOrtho_C
)
||
(
btype1
==
LibUtilities
::
eModified_B
)
||
(
btype1
==
LibUtilities
::
eModified_C
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 1 is for 2 or 3D expansions"
);
}
if
((
btype2
==
LibUtilities
::
eOrtho_B
)
||
(
btype2
==
LibUtilities
::
eOrtho_C
)
||
(
btype2
==
LibUtilities
::
eModified_B
)
||
(
btype2
==
LibUtilities
::
eModified_C
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 2 is for 2 or 3D expansions"
);
case
LibUtilitie
s
::
eQuadrilateral
:
if
((
btype1
==
LibUtilities
::
eOrtho_B
)
||
(
btype1
==
LibUtilities
::
eOrtho_C
)
||
(
btype1
==
LibUtilities
::
eModified_B
)
||
(
btype1
==
LibUtilities
::
eModified_C
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 1 is for 2 or 3D expansions"
);
}
if
((
btype2
==
LibUtilities
::
eOrtho_B
)
||
(
btype2
==
LibUtilities
::
eOrtho_C
)
||
(
btype2
==
LibUtilities
::
eModified_B
)
||
(
btype2
==
LibUtilities
::
eModified_C
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 2 is for 2 or 3D expansions"
);
}
break
;
default:
ASSERTL0
(
false
,
"Not a 2D expansion."
);
break
;
break
;
default:
ASSERTL0
(
false
,
"Not a 2D expansion."
);
break
;
}
order1
=
atoi
(
argv
[
4
]);
order2
=
atoi
(
argv
[
5
]);
nq1
=
atoi
(
argv
[
6
]);
...
...
@@ -161,7 +162,8 @@ int main(int argc, char *argv[])
if
(
btype2
!=
LibUtilities
::
eFourier
)
{
if
(
regionshape
==
StdRegions
::
eTriangle
)
{
if
(
regionshape
==
LibUtilities
::
eTriangle
)
{
Qtype2
=
LibUtilities
::
eGaussRadauMAlpha1Beta0
;
}
else
...
...
@@ -179,9 +181,9 @@ int main(int argc, char *argv[])
switch
(
regionshape
)
{
case
StdRegion
s
::
eTriangle
:
case
LibUtilitie
s
::
eTriangle
:
{
coords
[
0
]
=
atof
(
argv
[
8
]);
coords
[
1
]
=
atof
(
argv
[
9
]);
coords
[
2
]
=
atof
(
argv
[
10
]);
...
...
@@ -239,76 +241,76 @@ int main(int argc, char *argv[])
}
break
;
case
StdRegions
::
eQuadrilateral
:
{
// Gather coordinates
coords
[
0
]
=
atof
(
argv
[
8
]);
coords
[
1
]
=
atof
(
argv
[
9
]);
coords
[
2
]
=
atof
(
argv
[
10
]);
coords
[
3
]
=
atof
(
argv
[
11
]);
coords
[
4
]
=
atof
(
argv
[
12
]);
coords
[
5
]
=
atof
(
argv
[
13
]);
coords
[
6
]
=
atof
(
argv
[
14
]);
coords
[
7
]
=
atof
(
argv
[
15
]);
// Set up coordinates
const
int
zero
=
0
;
const
int
one
=
1
;
const
int
two
=
2
;
const
int
three
=
3
;
const
double
dZero
=
0.0
;
SpatialDomains
::
VertexComponentSharedPtr
verts
[
4
];
verts
[
0
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
zero
,
coords
[
0
],
coords
[
1
],
dZero
);
verts
[
1
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
one
,
coords
[
2
],
coords
[
3
],
dZero
);
verts
[
2
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
two
,
coords
[
4
],
coords
[
5
],
dZero
);
verts
[
3
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
three
,
coords
[
6
],
coords
[
7
],
dZero
);
// Set up Edges
SpatialDomains
::
SegGeomSharedPtr
edges
[
4
];
edges
[
0
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
zero
,
verts
[
0
],
verts
[
1
]);
edges
[
1
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
one
,
verts
[
1
],
verts
[
2
]);
edges
[
2
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
two
,
verts
[
2
],
verts
[
3
]);
edges
[
3
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
three
,
verts
[
3
],
verts
[
0
]);
StdRegions
::
Orientation
eorient
[
4
];
eorient
[
0
]
=
edgeDir
;
eorient
[
1
]
=
edgeDir
;
eorient
[
2
]
=
edgeDir
;
eorient
[
3
]
=
edgeDir
;
SpatialDomains
::
QuadGeomSharedPtr
geom
=
MemoryManager
<
SpatialDomains
::
QuadGeom
>::
AllocateSharedPtr
(
zero
,
verts
,
edges
,
eorient
);
geom
->
SetOwnData
();
const
LibUtilities
::
PointsKey
Pkey1
(
nq1
,
Qtype1
);
const
LibUtilities
::
PointsKey
Pkey2
(
nq2
,
Qtype2
);
const
LibUtilities
::
BasisKey
Bkey1
(
btype1
,
order1
,
Pkey1
);
const
LibUtilities
::
BasisKey
Bkey2
(
btype2
,
order2
,
Pkey2
);
E
=
new
LocalRegions
::
QuadExp
(
Bkey1
,
Bkey2
,
geom
);
//----------------------------------------------
// Define solution to be projected
E
->
GetCoords
(
x
,
y
);
for
(
i
=
0
;
i
<
nq1
*
nq2
;
++
i
)
case
LibUtilities
::
eQuadrilateral
:
{
sol
[
i
]
=
Quad_sol
(
x
[
i
],
y
[
i
],
order1
,
order2
,
btype1
,
btype2
);
// Gather coordinates
coords
[
0
]
=
atof
(
argv
[
8
]);
coords
[
1
]
=
atof
(
argv
[
9
]);
coords
[
2
]
=
atof
(
argv
[
10
]);
coords
[
3
]
=
atof
(
argv
[
11
]);
coords
[
4
]
=
atof
(
argv
[
12
]);
coords
[
5
]
=
atof
(
argv
[
13
]);
coords
[
6
]
=
atof
(
argv
[
14
]);
coords
[
7
]
=
atof
(
argv
[
15
]);
// Set up coordinates
const
int
zero
=
0
;
const
int
one
=
1
;
const
int
two
=
2
;
const
int
three
=
3
;
const
double
dZero
=
0.0
;
SpatialDomains
::
VertexComponentSharedPtr
verts
[
4
];
verts
[
0
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
zero
,
coords
[
0
],
coords
[
1
],
dZero
);
verts
[
1
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
one
,
coords
[
2
],
coords
[
3
],
dZero
);
verts
[
2
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
two
,
coords
[
4
],
coords
[
5
],
dZero
);
verts
[
3
]
=
MemoryManager
<
SpatialDomains
::
VertexComponent
>::
AllocateSharedPtr
(
two
,
three
,
coords
[
6
],
coords
[
7
],
dZero
);
// Set up Edges
SpatialDomains
::
SegGeomSharedPtr
edges
[
4
];
edges
[
0
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
zero
,
verts
[
0
],
verts
[
1
]);
edges
[
1
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
one
,
verts
[
1
],
verts
[
2
]);
edges
[
2
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
two
,
verts
[
2
],
verts
[
3
]);
edges
[
3
]
=
MemoryManager
<
SpatialDomains
::
SegGeom
>::
AllocateSharedPtr
(
three
,
verts
[
3
],
verts
[
0
]);
StdRegions
::
Orientation
eorient
[
4
];
eorient
[
0
]
=
edgeDir
;
eorient
[
1
]
=
edgeDir
;
eorient
[
2
]
=
edgeDir
;
eorient
[
3
]
=
edgeDir
;
SpatialDomains
::
QuadGeomSharedPtr
geom
=
MemoryManager
<
SpatialDomains
::
QuadGeom
>::
AllocateSharedPtr
(
zero
,
verts
,
edges
,
eorient
);
geom
->
SetOwnData
();
const
LibUtilities
::
PointsKey
Pkey1
(
nq1
,
Qtype1
);
const
LibUtilities
::
PointsKey
Pkey2
(
nq2
,
Qtype2
);
const
LibUtilities
::
BasisKey
Bkey1
(
btype1
,
order1
,
Pkey1
);
const
LibUtilities
::
BasisKey
Bkey2
(
btype2
,
order2
,
Pkey2
);
E
=
new
LocalRegions
::
QuadExp
(
Bkey1
,
Bkey2
,
geom
);
//----------------------------------------------
// Define solution to be projected
E
->
GetCoords
(
x
,
y
);
for
(
i
=
0
;
i
<
nq1
*
nq2
;
++
i
)
{
sol
[
i
]
=
Quad_sol
(
x
[
i
],
y
[
i
],
order1
,
order2
,
btype1
,
btype2
);
}
//---------------------------------------------
}
//---------------------------------------------
}
break
;
default:
ASSERTL0
(
false
,
"Not a 2D expansion."
);
break
;
default:
ASSERTL0
(
false
,
"Not a 2D expansion."
);
break
;
}
//--------------------------------------------
// Take the numerical derivative of the solution and add together in sol
E
->
PhysDeriv
(
sol
,
dx
,
dy
);
Vmath
::
Vadd
(
nq1
*
nq2
,
dx
,
1
,
dy
,
1
,
sol
,
1
);
//---------------------------------------------
//---------------------------------------------
// Project onto Expansion
E
->
FwdTrans
(
sol
,
E
->
UpdateCoeffs
());
...
...
@@ -323,7 +325,7 @@ int main(int argc, char *argv[])
// Define exact solution of differential
switch
(
regionshape
)
{
case
StdRegion
s
::
eTriangle
:
case
LibUtilitie
s
::
eTriangle
:
{
//----------------------------------------------
// Define solution to be differentiated
...
...
@@ -334,7 +336,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
}
break
;
case
StdRegion
s
::
eQuadrilateral
:
case
LibUtilitie
s
::
eQuadrilateral
:
{
for
(
i
=
0
;
i
<
nq1
*
nq2
;
++
i
)
{
...
...
@@ -343,11 +345,11 @@ int main(int argc, char *argv[])
}
//---------------------------------------------
break
;
default:
ASSERTL0
(
false
,
"Not a 2D expansion."
);
break
;
default:
ASSERTL0
(
false
,
"Not a 2D expansion."
);
break
;
}
//--------------------------------------------
// Write solution
ofstream
outfile
(
"ProjectFile2D.dat"
);
...
...
library/Demos/LocalRegions/LocProject_Diff3D.cpp
View file @
612a6e36
...
...
@@ -73,7 +73,7 @@ int main(int argc, char *argv[]){
LibUtilities
::
PointsType
Qtype1
,
Qtype2
,
Qtype3
;
LibUtilities
::
BasisType
btype1
,
btype2
,
btype3
;
StdRegions
::
Expansion
Type
regionshape
;
LibUtilities
::
Shape
Type
regionshape
;
StdRegions
::
StdExpansion
*
E
;
Array
<
OneD
,
NekDouble
>
x
,
y
,
z
,
sol
,
dx
,
dy
,
dz
;
...
...
@@ -85,9 +85,9 @@ int main(int argc, char *argv[]){
fprintf
(
stderr
,
"Where RegionShape is an integer value which "
"dictates the region shape:
\n
"
);
fprintf
(
stderr
,
"
\t
Tetrahedron =
4
\n
"
);
fprintf
(
stderr
,
"
\t
Prism =
6
\n
"
);
fprintf
(
stderr
,
"
\t
Hexahedron =
7
\n
"
);
fprintf
(
stderr
,
"
\t
Tetrahedron =
5
\n
"
);
fprintf
(
stderr
,
"
\t
Prism =
7
\n
"
);
fprintf
(
stderr
,
"
\t
Hexahedron =
8
\n
"
);
fprintf
(
stderr
,
"Where type is an integer value which "
"dictates the basis as:
\n
"
);
...
...
@@ -106,12 +106,12 @@ int main(int argc, char *argv[]){
exit
(
1
);
}
regionshape
=
(
Expansion
Type
)
atoi
(
argv
[
1
]);
regionshape
=
(
Shape
Type
)
atoi
(
argv
[
1
]);
// Check to see if 3D region
if
(
regionshape
!=
StdRegion
s
::
eTetrahedron
&&
regionshape
!=
StdRegion
s
::
ePrism
&&
regionshape
!=
StdRegion
s
::
eHexahedron
)
if
(
regionshape
!=
LibUtilitie
s
::
eTetrahedron
&&
regionshape
!=
LibUtilitie
s
::
ePrism
&&
regionshape
!=
LibUtilitie
s
::
eHexahedron
)
{
NEKERROR
(
ErrorUtil
::
efatal
,
"This shape is not a 3D region"
);
}
...
...
@@ -127,74 +127,74 @@ int main(int argc, char *argv[]){
// Check to see that correct Expansions are used
switch
(
regionshape
)
{
case
StdRegion
s
::
eTetrahedron
:
if
((
btype1
==
eOrtho_B
)
||
(
btype1
==
eOrtho_C
)
||
(
btype1
==
eModified_B
)
||
(
btype1
==
eModified_C
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 1 cannot be of type Ortho_B, "
"Ortho_C, Modified_B or Modified_C"
);
}
if
((
btype2
==
eOrtho_A
)
||
(
btype2
==
eOrtho_C
)
||
(
btype2
==
eModified_A
)
||
(
btype2
==
eModified_C
))
{
NEKERROR
(
ErrorUtil
::
efatal
,
"Basis 2 cannot be of type Ortho_A, "
"Ortho_C, Modified_A or Modified_C"
);
}
if
((
btype3
==
eOrtho_A
)
||
(
btype3
==
eOrtho_B
)
||
(
btype3
==
eModified_A
)
||
(
bt