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
Julia
Nektar
Commits
27a6ee0c
Commit
27a6ee0c
authored
Jul 06, 2017
by
Julia Docampo Sanchez
Browse files
compiling version: TP FILTERS and 2D LINE FILTER
parent
e77133a4
Changes
157
Expand all
Hide whitespace changes
Inline
Side-by-side
library/CMakeLists.txt
View file @
27a6ee0c
SET
(
LibrarySubDirs FieldUtils GlobalMapping LibUtilities LocalRegions
Collections MultiRegions SpatialDomains StdRegions
L
SIAC LSIAC
2D
)
Collections MultiRegions SpatialDomains StdRegions
TP
SIAC LSIAC
)
SET
(
UnitTestSubDirs UnitTests
)
SET
(
DemoSubDirs Demos
)
SET
(
TimingsSubDirs Timings
)
...
...
library/LSIAC/BSplines.h
View file @
27a6ee0c
...
...
@@ -3,7 +3,7 @@
#include <iostream>
#include <vector>
#include "
NektarBaseClass
.h"
#include "
LSIACPostProcessor
.h"
using
namespace
std
;
/// This class is the base class for all the B-Splines.
...
...
@@ -17,7 +17,7 @@ namespace LSIAC
{
class
BSplines
:
public
NektarBaseClass
{
class
BSplines
:
public
LSIACPostProcessor
{
protected:
public:
};
...
...
library/LSIAC/CMakeLists.txt
View file @
27a6ee0c
SET
(
LSIACHeaders
NektarBaseClass.h
BSplines.h
FilterType.h
LSIACPostProcessor.h
SIACenumerators.h
GeneralBSplines.h
CentralBSplines.h
SIAC
Filter
.h
Symmetric
SIAC
.h
OneSided
SIAC
.h
SIAC
Kernel
.h
Symmetric
Kernel
.h
OneSided
Kernel
.h
HandleMesh.h
HandleNekMesh.h
HandleNekMesh1D.h
HandleNekMesh2D.h
HandleNekMesh3D.h
Smoothie.h
SmoothieSIAC.h
SmoothieSIAC1D.h
SmoothieSIAC2D.h
SmoothieSIAC3D.h
element_centre_size.h
SIACpostprocessor.h
SIACFilter.h
SIACFilter1D.h
SIACFilter2D.h
SIACFilter3D.h
)
SET
(
LSIACSources
NektarBaseClass.cpp
fieldDerivatives.cpp
LSIACPostProcessor.cpp
GeneralBSplines.cpp
CentralBSplines.cpp
SIAC
Filter
.cpp
Symmetric
SIAC
.cpp
OneSided
SIAC
.cpp
SIAC
Kernel
.cpp
Symmetric
Kernel
.cpp
OneSided
Kernel
.cpp
HandleNekMesh.cpp
HandleNekMesh1D.cpp
HandleNekMesh2D.cpp
HandleNekMesh3D.cpp
SmoothieSIAC.cpp
SmoothieSIAC1D.cpp
SmoothieSIAC2D.cpp
SmoothieSIAC3D.cpp
SIACFilter.cpp
SIACFilter1D.cpp
SIACFilter2D.cpp
SIACFilter3D.cpp
ppio.cpp
)
ADD_NEKTAR_LIBRARY
(
LSIAC lib
${
NEKTAR_LIBRARY_TYPE
}
${
LSIACSources
}
${
LSIACHeaders
}
)
...
...
library/LSIAC/CMakeLists.txt.user
0 → 100644
View file @
27a6ee0c
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE QtCreatorProject>
<!-- Written by QtCreator 3.0.1, 2017-07-05T11:26:18. -->
<qtcreator>
<data>
<variable>
ProjectExplorer.Project.ActiveTarget
</variable>
<value
type=
"int"
>
-1
</value>
</data>
<data>
<variable>
ProjectExplorer.Project.EditorSettings
</variable>
<valuemap
type=
"QVariantMap"
>
<value
type=
"bool"
key=
"EditorConfiguration.AutoIndent"
>
true
</value>
<value
type=
"bool"
key=
"EditorConfiguration.AutoSpacesForTabs"
>
false
</value>
<value
type=
"bool"
key=
"EditorConfiguration.CamelCaseNavigation"
>
true
</value>
<valuemap
type=
"QVariantMap"
key=
"EditorConfiguration.CodeStyle.0"
>
<value
type=
"QString"
key=
"language"
>
Cpp
</value>
<valuemap
type=
"QVariantMap"
key=
"value"
>
<value
type=
"QByteArray"
key=
"CurrentPreferences"
>
CppGlobal
</value>
</valuemap>
</valuemap>
<valuemap
type=
"QVariantMap"
key=
"EditorConfiguration.CodeStyle.1"
>
<value
type=
"QString"
key=
"language"
>
QmlJS
</value>
<valuemap
type=
"QVariantMap"
key=
"value"
>
<value
type=
"QByteArray"
key=
"CurrentPreferences"
>
QmlJSGlobal
</value>
</valuemap>
</valuemap>
<value
type=
"int"
key=
"EditorConfiguration.CodeStyle.Count"
>
2
</value>
<value
type=
"QByteArray"
key=
"EditorConfiguration.Codec"
>
UTF-8
</value>
<value
type=
"bool"
key=
"EditorConfiguration.ConstrainTooltips"
>
false
</value>
<value
type=
"int"
key=
"EditorConfiguration.IndentSize"
>
4
</value>
<value
type=
"bool"
key=
"EditorConfiguration.KeyboardTooltips"
>
false
</value>
<value
type=
"bool"
key=
"EditorConfiguration.MouseNavigation"
>
true
</value>
<value
type=
"int"
key=
"EditorConfiguration.PaddingMode"
>
1
</value>
<value
type=
"bool"
key=
"EditorConfiguration.ScrollWheelZooming"
>
true
</value>
<value
type=
"int"
key=
"EditorConfiguration.SmartBackspaceBehavior"
>
1
</value>
<value
type=
"bool"
key=
"EditorConfiguration.SpacesForTabs"
>
true
</value>
<value
type=
"int"
key=
"EditorConfiguration.TabKeyBehavior"
>
2
</value>
<value
type=
"int"
key=
"EditorConfiguration.TabSize"
>
8
</value>
<value
type=
"bool"
key=
"EditorConfiguration.UseGlobal"
>
true
</value>
<value
type=
"int"
key=
"EditorConfiguration.Utf8BomBehavior"
>
1
</value>
<value
type=
"bool"
key=
"EditorConfiguration.addFinalNewLine"
>
true
</value>
<value
type=
"bool"
key=
"EditorConfiguration.cleanIndentation"
>
true
</value>
<value
type=
"bool"
key=
"EditorConfiguration.cleanWhitespace"
>
true
</value>
<value
type=
"bool"
key=
"EditorConfiguration.inEntireDocument"
>
false
</value>
</valuemap>
</data>
<data>
<variable>
ProjectExplorer.Project.PluginSettings
</variable>
<valuemap
type=
"QVariantMap"
/>
</data>
<data>
<variable>
ProjectExplorer.Project.TargetCount
</variable>
<value
type=
"int"
>
0
</value>
</data>
<data>
<variable>
ProjectExplorer.Project.Updater.EnvironmentId
</variable>
<value
type=
"QByteArray"
>
{a578d475-18b0-40f9-843f-5eee755f73a3}
</value>
</data>
<data>
<variable>
ProjectExplorer.Project.Updater.FileVersion
</variable>
<value
type=
"int"
>
15
</value>
</data>
</qtcreator>
library/LSIAC/FilterType.h
deleted
100644 → 0
View file @
e77133a4
#ifndef FILTERTYPE_H
#define FILTERTYPE_H
namespace
Nektar
{
namespace
LSIAC
{
namespace
SIACUtilities
{
//! This enum specifies the kind of filter to be used by the object.
/*! The user can specify any type of smoothing or derivative filter.
* All the properties of memeber depends on this enum. This enum cannot be changed after the object has been initialized.
* A new enum needs to be added to the list to create any new type of filter.
*/
enum
FilterType
{
eNONE
,
/// 0
ePURE
,
ePURE_2kp1
,
eSTD
,
/// 1
eSTD_1SIDED_XLI
,
/// 2
eSTD_2kp1
,
/// 3
eSTD_2kp1_1SIDED_XLI
,
/// 4
eSTD_2kp1_1SIDED_4kp1
,
/// 5
eSTD_2kp1_1SIDED_SRV
,
/// 6
eDER
,
/// 7
eDER_1SIDED
,
///8
eDER_1SIDED_XLI
,
/// 9
eDER_2kp1
,
/// 10
eDER_2kp1_1SIDED
,
/// 11
eDER_2kp1_1SIDED_XLI
,
/// 12
eDER_2kp1_1SIDED_4kp1
,
/// 13
eDER_2kp1_1SIDED_SRV
,
/// 14
eDER_LOWER_ORDER
,
eDER_2kp1_LOWER_ORDER
};
inline
void
write_filter_name
(
FilterType
filter
,
string
&
name
)
{
switch
(
filter
)
{
case
eNONE
:
name
=
"NONE"
;
return
;
case
ePURE
:
name
=
"PURE"
;
return
;
case
ePURE_2kp1
:
name
=
"PURE_2kp1"
;
return
;
case
eSTD
:
name
=
"STD"
;
return
;
case
eSTD_1SIDED_XLI
:
name
=
"STD_1SIDED_XLI"
;
return
;
case
eSTD_2kp1
:
name
=
"STD_2kp1"
;
return
;
case
eSTD_2kp1_1SIDED_XLI
:
name
=
"STD_2kp1_1SIDED_XLI"
;
return
;
case
eSTD_2kp1_1SIDED_4kp1
:
name
=
"STD_2kp1_1SIDED_4kp1"
;
return
;
case
eSTD_2kp1_1SIDED_SRV
:
name
=
"STD_2kp1_1SIDED_SRV"
;
return
;
case
eDER
:
name
=
"DER"
;
return
;
case
eDER_1SIDED
:
name
=
"DER_1SIDED"
;
return
;
case
eDER_1SIDED_XLI
:
name
=
"DER_1SIDED_XLI"
;
return
;
case
eDER_2kp1
:
name
=
"DER_2kp1"
;
return
;
case
eDER_2kp1_1SIDED
:
name
=
"DER_2kp1_1SIDED"
;
return
;
case
eDER_2kp1_1SIDED_XLI
:
name
=
"DER_2kp1_1SIDED_XLI"
;
return
;
case
eDER_2kp1_1SIDED_4kp1
:
name
=
"DER_2kp1_1SIDED_4kp1"
;
return
;
case
eDER_2kp1_1SIDED_SRV
:
name
=
"DER_2kp1_1SIDED_SRV"
;
return
;
case
eDER_LOWER_ORDER
:
name
=
"DER_LOWER_ORDER"
;
return
;
case
eDER_2kp1_LOWER_ORDER
:
name
=
"DER_2kp1_LOWER_ORDER"
;
return
;
}
}
struct
SMOOTHIE_Status
{
int
count
;
int
cancelled
;
int
SOURCE
;
int
TAG
;
int
ERROR
;
};
}
}
}
#endif
library/LSIAC/HandleMesh.h
View file @
27a6ee0c
...
...
@@ -4,7 +4,7 @@
#include <iostream>
#include <vector>
#include <string>
#include "
NektarBaseClass
.h"
#include "
LSIACPostProcessor
.h"
using
namespace
std
;
...
...
@@ -19,20 +19,10 @@ namespace Nektar
namespace
LSIAC
{
class
HandleMesh
:
public
NektarBaseClass
{
class
HandleMesh
:
public
LSIACPostProcessor
{
public:
protected:
};
}
...
...
library/LSIAC/HandleNekMesh.cpp
View file @
27a6ee0c
This diff is collapsed.
Click to expand it.
library/LSIAC/HandleNekMesh.h
View file @
27a6ee0c
This diff is collapsed.
Click to expand it.
library/LSIAC/HandleNekMesh1D.h
View file @
27a6ee0c
...
...
@@ -20,11 +20,6 @@ namespace Nektar
{
m_useRTree
=
false
;
// cout << "into HandleNekMesh1D constructor" << endl;
// m_expansions.push_back(MemoryManager<MultiRegions::ContField1D>::
// AllocateSharedPtr(m_session,m_graph, m_session->GetVariable(0)) );
// cout << "Expansion Type: " << m_expansions[0]->GetExpType() << endl;
}
protected:
...
...
library/LSIAC/HandleNekMesh2D.cpp
View file @
27a6ee0c
...
...
@@ -6,146 +6,154 @@
#include <iomanip> // std::setprecision
#include <limits>
namespace
Nektar
{
namespace
LSIAC
{
bool
HandleNekMesh2D
::
v_LoadMesh
(
string
var
)
{
SpatialDomains
::
ExpansionMap
expansions
=
m_graph
->
GetExpansions
();
m_expansions
.
push_back
(
MemoryManager
<
MultiRegions
::
DisContField2D
>
::
AllocateSharedPtr
(
m_session
,
m_graph
,
var
));
if
(
m_useRTree
)
LoadExpListIntoRTree
();
return
true
;
}
bool
HandleNekMesh2D
::
v_LoadData
(
string
filename
,
vector
<
string
>
&
variables
)
{
SpatialDomains
::
ExpansionMap
expansions
=
m_graph
->
GetExpansions
();
for
(
int
i
=
0
;
i
<
variables
.
size
();
i
++
)
{
m_expansions
.
push_back
(
MemoryManager
<
MultiRegions
::
DisContField2D
>
::
AllocateSharedPtr
(
m_session
,
m_graph
,
variables
[
i
]
));
}
std
::
vector
<
LibUtilities
::
FieldDefinitionsSharedPtr
>
rFieldDef
;
std
::
vector
<
std
::
vector
<
NekDouble
>
>
rFieldData
;
Array
<
OneD
,
int
>
ElementGIDs
(
expansions
.
size
());
SpatialDomains
::
ExpansionMap
::
const_iterator
expIt
;
int
i
=
0
;
for
(
expIt
=
expansions
.
begin
();
expIt
!=
expansions
.
end
();
++
expIt
)
{
ElementGIDs
[
i
++
]
=
expIt
->
second
->
m_geomShPtr
->
GetGlobalID
();
}
Import
(
filename
,
rFieldDef
,
rFieldData
,
LibUtilities
::
NullFieldMetaDataMap
,
ElementGIDs
);
for
(
int
i
=
0
;
i
<
rFieldDef
.
size
()
;
i
++
)
{
for
(
int
e
=
0
;
e
<
variables
.
size
();
e
++
)
{
m_expansions
[
e
]
->
ExtractDataToCoeffs
(
rFieldDef
[
i
],
rFieldData
[
i
],
variables
[
e
],
m_expansions
[
e
]
->
UpdateCoeffs
());
}
}
for
(
int
i
=
0
;
i
<
m_expansions
.
size
();
i
++
)
{
m_expansions
[
i
]
->
BwdTrans
(
m_expansions
[
i
]
->
GetCoeffs
(),
m_expansions
[
i
]
->
UpdatePhys
());
}
return
true
;
}
bool
HandleNekMesh2D
::
v_GetBreakPts
(
const
Array
<
OneD
,
NekDouble
>
point_offset
,
Array
<
OneD
,
NekDouble
>
&
direction
,
const
NekDouble
tmin
,
const
NekDouble
tmax
,
vector
<
NekDouble
>
&
tPos
)
{
return
IntersectWithEdges
(
tPos
,
point_offset
,
direction
,
tmin
,
tmax
);
}
bool
HandleNekMesh2D
::
v_intersect
(
Array
<
OneD
,
NekDouble
>
&
p1
,
Array
<
OneD
,
NekDouble
>
&
p2
,
Array
<
OneD
,
NekDouble
>
&
q1
,
Array
<
OneD
,
NekDouble
>
&
q2
,
Array
<
OneD
,
NekDouble
>&
i1
,
vector
<
NekDouble
>&
i2
)
{
Array
<
OneD
,
NekDouble
>
threeDp1
(
threeD
,
0.
),
threeDp2
(
threeD
,
0.
),
threeDq1
(
threeD
,
0.
),
threeDq2
(
threeD
,
0.
);
i1
=
Array
<
OneD
,
NekDouble
>
(
threeD
,
0.0
);
for
(
int
d
=
0
;
d
<
p1
.
num_elements
();
++
d
)
{
threeDp1
[
d
]
=
p1
[
d
];
threeDp2
[
d
]
=
p2
[
d
];
threeDq1
[
d
]
=
q1
[
d
];
threeDq2
[
d
]
=
q2
[
d
];
}
return
(
threeD_intersect
(
threeDp1
,
threeDp2
,
threeDq1
,
threeDq2
,
i1
,
i2
));
}
{
namespace
LSIAC
{
bool
HandleNekMesh2D
::
v_LoadMesh
(
string
var
)
{
SpatialDomains
::
ExpansionMap
expansions
=
m_graph
->
GetExpansions
();
m_expansions
.
push_back
(
MemoryManager
<
MultiRegions
::
DisContField2D
>
::
AllocateSharedPtr
(
m_session
,
m_graph
,
var
));
if
(
m_useRTree
)
LoadExpListIntoRTree
();
return
true
;
}
bool
HandleNekMesh2D
::
v_LoadData
(
string
filename
,
vector
<
string
>
&
variables
)
{
SpatialDomains
::
ExpansionMap
expansions
=
m_graph
->
GetExpansions
();
for
(
int
i
=
0
;
i
<
variables
.
size
();
i
++
)
{
m_expansions
.
push_back
(
MemoryManager
<
MultiRegions
::
DisContField2D
>
::
AllocateSharedPtr
(
m_session
,
m_graph
,
variables
[
i
]
));
}
if
(
m_useRTree
)
LoadExpListIntoRTree
();
std
::
vector
<
LibUtilities
::
FieldDefinitionsSharedPtr
>
rFieldDef
;
std
::
vector
<
std
::
vector
<
NekDouble
>
>
rFieldData
;
Array
<
OneD
,
int
>
ElementGIDs
(
expansions
.
size
());
SpatialDomains
::
ExpansionMap
::
const_iterator
expIt
;
int
i
=
0
;
for
(
expIt
=
expansions
.
begin
();
expIt
!=
expansions
.
end
();
++
expIt
)
{
ElementGIDs
[
i
++
]
=
expIt
->
second
->
m_geomShPtr
->
GetGlobalID
();
}
Import
(
filename
,
rFieldDef
,
rFieldData
,
LibUtilities
::
NullFieldMetaDataMap
,
ElementGIDs
);
for
(
int
i
=
0
;
i
<
rFieldDef
.
size
()
;
i
++
)
{
for
(
int
e
=
0
;
e
<
variables
.
size
();
e
++
)
{
m_expansions
[
e
]
->
ExtractDataToCoeffs
(
rFieldDef
[
i
],
rFieldData
[
i
],
variables
[
e
],
m_expansions
[
e
]
->
UpdateCoeffs
());
}
}
for
(
int
i
=
0
;
i
<
m_expansions
.
size
();
i
++
)
{
m_expansions
[
i
]
->
BwdTrans
(
m_expansions
[
i
]
->
GetCoeffs
(),
m_expansions
[
i
]
->
UpdatePhys
());
}
return
true
;
}
// Mostly funtions using RTree from here.
int
HandleNekMesh2D
::
v_GetExpansionIndexUsingRTree
(
const
Array
<
OneD
,
NekDouble
>
&
point
)
const
{
std
::
vector
<
RValue
>
result_s
;
m_rtree
.
query
(
Boostgi
::
nearest
(
RPoint
(
point
[
0
],
point
[
1
]),
1
),
std
::
back_inserter
(
result_s
));
assert
(
result_s
.
size
()
==
1
&&
"no element found. Something is cleraly wrong"
);
// Can write better code using within.
RValue
v
=
result_s
[
0
];
return
boost
::
get
<
1
>
(
v
);
// <2> is global ID.
}
void
HandleNekMesh2D
::
GetBoundingOfElement
(
SpatialDomains
::
GeometrySharedPtr
elGeom
,
RBox
&
b
)
{
int
dim
=
twoD
;
Array
<
OneD
,
NekDouble
>
vertex_coord
(
dim
),
MINP
(
dim
,
std
::
numeric_limits
<
NekDouble
>::
max
()),
MAXP
(
dim
,
std
::
numeric_limits
<
NekDouble
>::
min
());
int
numPts
=
elGeom
->
GetNumVerts
();
for
(
int
v
=
0
;
v
<
numPts
;
v
++
)
{
elGeom
->
GetVertex
(
v
)
->
GetCoords
(
vertex_coord
);
for
(
int
d
=
0
;
d
<
dim
;
++
d
)
{
MINP
[
d
]
=
std
::
min
(
MINP
[
d
],
vertex_coord
[
d
]);
MAXP
[
d
]
=
std
::
max
(
MAXP
[
d
],
vertex_coord
[
d
]);
}
}
Boostg
::
set
<
Boostg
::
min_corner
,
0
>
(
b
,
MINP
[
0
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
min_corner
,
1
>
(
b
,
MINP
[
1
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
0
>
(
b
,
MAXP
[
0
]
+
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
1
>
(
b
,
MAXP
[
1
]
+
TOLERENCE_MESH_COMP
);
return
;
}
void
HandleNekMesh2D
::
LoadExpListIntoRTree
()
{
// check if expansions are already loaded.
assert
(
m_expansions
.
size
()
>
0
&&
"should have loaded atleast one expansion"
);
MultiRegions
::
ExpListSharedPtr
expList
=
m_expansions
[
0
];
int
expSize
=
expList
->
GetExpSize
();
for
(
int
e
=
0
;
e
<
expSize
;
++
e
)
{
SpatialDomains
::
GeometrySharedPtr
geom
=
expList
->
GetExp
(
e
)
->
GetGeom
();
int
gid
=
geom
->
GetGlobalID
();
RBox
b
;
GetBoundingOfElement
(
geom
,
b
);
m_rtree
.
insert
(
boost
::
make_tuple
(
b
,
e
,
gid
)
);
}
m_useRTree
=
true
;
}
void
HandleNekMesh2D
::
IntersectWithBoxUsingRTree
(
Array
<
OneD
,
NekDouble
>
&
minCorner
,
Array
<
OneD
,
NekDouble
>
maxCorner
,
vector
<
int
>
&
elIds
,
vector
<
int
>
&
glIds
)
{
RBox
b
;
Boostg
::
set
<
Boostg
::
min_corner
,
0
>
(
b
,
minCorner
[
0
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
min_corner
,
1
>
(
b
,
minCorner
[
1
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
0
>
(
b
,
maxCorner
[
0
]
+
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
1
>
(
b
,
maxCorner
[
1
]
+
TOLERENCE_MESH_COMP
);
//rtree.query()
std
::
vector
<
RValue
>
result_s
;
m_rtree
.
query
(
Boostgi
::
intersects
(
b
),
std
::
back_inserter
(
result_s
)
);
elIds
.
clear
();
glIds
.
clear
();
//assert( false && "Not done developing yet :( ");
BOOST_FOREACH
(
RValue
const
&
v
,
result_s
)
{
int
eId
=
boost
::
get
<
1
>
(
v
);
int
gId
=
boost
::
get
<
2
>
(
v
);
elIds
.
push_back
(
eId
);
glIds
.
push_back
(
gId
);
}
}
}
}
bool
HandleNekMesh2D
::
v_GetBreakPts
(
const
Array
<
OneD
,
NekDouble
>
point_offset
,
Array
<
OneD
,
NekDouble
>
&
direction
,
const
NekDouble
tmin
,
const
NekDouble
tmax
,
vector
<
NekDouble
>
&
tPos
)
{
return
IntersectWithEdges
(
tPos
,
point_offset
,
direction
,
tmin
,
tmax
);
}
bool
HandleNekMesh2D
::
v_intersect
(
Array
<
OneD
,
NekDouble
>
&
p1
,
Array
<
OneD
,
NekDouble
>
&
p2
,
Array
<
OneD
,
NekDouble
>
&
q1
,
Array
<
OneD
,
NekDouble
>
&
q2
,
Array
<
OneD
,
NekDouble
>&
i1
,
vector
<
NekDouble
>&
i2
)
{
int
three
=
3
;
Array
<
OneD
,
NekDouble
>
threeDp1
(
three
,
0.
),
threeDp2
(
three
,
0.
),
threeDq1
(
three
,
0.
),
threeDq2
(
three
,
0.
);
i1
=
Array
<
OneD
,
NekDouble
>
(
three
,
0.0
);
for
(
int
d
=
0
;
d
<
p1
.
num_elements
();
++
d
)
{
threeDp1
[
d
]
=
p1
[
d
];
threeDp2
[
d
]
=
p2
[
d
];
threeDq1
[
d
]
=
q1
[
d
];
threeDq2
[
d
]
=
q2
[
d
];
}
return
(
threeD_intersect
(
threeDp1
,
threeDp2
,
threeDq1
,
threeDq2
,
i1
,
i2
));
}
// Mostly funtions using RTree from here.
int
HandleNekMesh2D
::
v_GetExpansionIndexUsingRTree
(
const
Array
<
OneD
,
NekDouble
>
&
point
)
const
{
std
::
vector
<
RValue
>
result_s
;
m_rtree
.
query
(
Boostgi
::
nearest
(
RPoint
(
point
[
0
],
point
[
1
]),
1
),
std
::
back_inserter
(
result_s
));
assert
(
result_s
.
size
()
==
1
&&
"no element found. Something is cleraly wrong"
);
// Can write better code using within.
RValue
v
=
result_s
[
0
];
return
boost
::
get
<
1
>
(
v
);
// <2> is global ID.
}
void
HandleNekMesh2D
::
GetBoundingOfElement
(
SpatialDomains
::
GeometrySharedPtr
elGeom
,
RBox
&
b
)
{
int
dim
=
2
;
Array
<
OneD
,
NekDouble
>
vertex_coord
(
dim
),
MINP
(
dim
,
std
::
numeric_limits
<
NekDouble
>::
max
()),
MAXP
(
dim
,
std
::
numeric_limits
<
NekDouble
>::
min
());
int
numPts
=
elGeom
->
GetNumVerts
();
for
(
int
v
=
0
;
v
<
numPts
;
v
++
)
{
elGeom
->
GetVertex
(
v
)
->
GetCoords
(
vertex_coord
);
for
(
int
d
=
0
;
d
<
dim
;
++
d
)
{
MINP
[
d
]
=
std
::
min
(
MINP
[
d
],
vertex_coord
[
d
]);
MAXP
[
d
]
=
std
::
max
(
MAXP
[
d
],
vertex_coord
[
d
]);
}
}
Boostg
::
set
<
Boostg
::
min_corner
,
0
>
(
b
,
MINP
[
0
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
min_corner
,
1
>
(
b
,
MINP
[
1
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
0
>
(
b
,
MAXP
[
0
]
+
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
1
>
(
b
,
MAXP
[
1
]
+
TOLERENCE_MESH_COMP
);
return
;
}
void
HandleNekMesh2D
::
LoadExpListIntoRTree
()
{
// check if expansions are already loaded.
assert
(
m_expansions
.
size
()
>
0
&&
"should have loaded atleast one expansion"
);
MultiRegions
::
ExpListSharedPtr
expList
=
m_expansions
[
0
];
int
expSize
=
expList
->
GetExpSize
();
for
(
int
e
=
0
;
e
<
expSize
;
++
e
)
{
SpatialDomains
::
GeometrySharedPtr
geom
=
expList
->
GetExp
(
e
)
->
GetGeom
();
int
gid
=
geom
->
GetGlobalID
();
RBox
b
;
GetBoundingOfElement
(
geom
,
b
);
m_rtree
.
insert
(
boost
::
make_tuple
(
b
,
e
,
gid
)
);
}
m_useRTree
=
true
;
}
void
HandleNekMesh2D
::
IntersectWithBoxUsingRTree
(
Array
<
OneD
,
NekDouble
>
&
minCorner
,
Array
<
OneD
,
NekDouble
>
maxCorner
,
vector
<
int
>
&
elIds
,
vector
<
int
>
&
glIds
)
{
RBox
b
;
Boostg
::
set
<
Boostg
::
min_corner
,
0
>
(
b
,
minCorner
[
0
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
min_corner
,
1
>
(
b
,
minCorner
[
1
]
-
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
0
>
(
b
,
maxCorner
[
0
]
+
TOLERENCE_MESH_COMP
);
Boostg
::
set
<
Boostg
::
max_corner
,
1
>
(
b
,
maxCorner
[
1
]
+
TOLERENCE_MESH_COMP
);
//rtree.query()
std
::
vector
<
RValue
>
result_s
;
m_rtree
.
query
(
Boostgi
::
intersects
(
b
),
std
::
back_inserter
(
result_s
)
);
elIds
.
clear
();
glIds
.
clear
();
//assert( false && "Not done developing yet :( ");
BOOST_FOREACH
(
RValue
const
&
v
,
result_s
)
{