48 namespace LibUtilities
55 bool isVertex(
int x,
int y,
int z,
int npts){
56 return (x==0 && y==0 && z==0) || (x==(npts-1) && y==0 && z==0) || (x==0 && y==(npts-1) && z==0) || (x==0 && y==0 && z==(npts-1));
59 bool isEdge_01(
int x,
int y,
int z,
int npts){
63 bool isEdge_12(
int x,
int y,
int z,
int npts){
64 return z==0 && x + y == npts -1;
67 bool isEdge_20(
int x,
int y,
int z,
int npts){
71 bool isEdge_03(
int x,
int y,
int z,
int npts){
75 bool isEdge_13(
int x,
int y,
int z,
int npts){
76 return y==0 && x + z == npts -1;
79 bool isEdge_23(
int x,
int y,
int z,
int npts){
80 return x==0 && y + z == npts -1;
83 bool isEdge(
int x,
int y,
int z,
int npts){
84 return isEdge_01(x,y,z,npts)||isEdge_12(x,y,z,npts)||isEdge_20(x,y,z,npts)
85 ||isEdge_03(x,y,z,npts)||isEdge_13(x,y,z,npts)||isEdge_23(x,y,z,npts);
88 bool isFace_012(
int x,
int y,
int z,
int npts){
92 bool isFace_013(
int x,
int y,
int z,
int npts){
96 bool isFace_123(
int x,
int y,
int z,
int npts){
100 bool isFace_203(
int x,
int y,
int z,
int npts){
104 bool isFace(
int x,
int y,
int z,
int npts){
105 return isFace_012(x, y, z, npts) || isFace_013(x, y, z, npts)
106 || isFace_123(x, y, z, npts) || isFace_203(x, y, z, npts);
119 for(
unsigned int z=0, index=0; z<
npts; ++z){
120 for(
int y=0; y<npts-z; ++y){
121 for(
int x=0; x<npts-z-y; ++x, ++index){
141 vector<int> iEdge_01;
142 vector<int> iEdge_12;
143 vector<int> iEdge_20;
144 vector<int> iEdge_03;
145 vector<int> iEdge_13;
146 vector<int> iEdge_23;
147 vector<int> iFace_012;
148 vector<int> iFace_013;
149 vector<int> iFace_123;
150 vector<int> iFace_203;
151 vector<int> interiorVolumePoints;
155 for(
int z=0, index=0; z<
npts; ++z){
156 for(
int y=0; y<npts-z; ++y){
157 for(
int x=0; x<npts-z-y; ++x, ++index){
159 if( isVertex(x,y,z,npts) ){
161 vertex.push_back(index);
163 }
else if( isEdge(x,y,z,npts) ){
165 if( isEdge_01(x,y,z,npts) ){
167 iEdge_01.push_back(index);
169 }
else if( isEdge_12(x,y,z,npts) ){
171 iEdge_12.push_back(index);
173 }
else if( isEdge_20(x,y,z,npts) ){
175 iEdge_20.insert(iEdge_20.begin(), index);
177 }
else if( isEdge_03(x,y,z,npts) ){
179 iEdge_03.push_back(index);
181 }
else if( isEdge_13(x,y,z,npts) ){
183 iEdge_13.push_back(index);
185 }
else if( isEdge_23(x,y,z,npts) ){
187 iEdge_23.push_back(index);
191 }
else if( isFace(x,y,z,npts) ) {
193 if( isFace_012(x,y,z,npts) ){
195 iFace_012.push_back(index);
197 }
else if( isFace_013(x,y,z,npts) ){
199 iFace_013.push_back(index);
201 }
else if( isFace_123(x,y,z,npts) ){
203 iFace_123.push_back(index);
205 }
else if( isFace_203(x,y,z,npts) ){
207 iFace_203.push_back(index);
212 interiorVolumePoints.push_back(index);
220 for(
unsigned int n=0; n<vertex.size(); ++n){
222 map.push_back(vertex[n]);
225 for(
unsigned int n=0; n<iEdge_01.size(); ++n){
227 map.push_back(iEdge_01[n]);
230 for(
unsigned int n=0; n<iEdge_12.size(); ++n){
232 map.push_back(iEdge_12[n]);
235 for(
unsigned int n=0; n<iEdge_20.size(); ++n){
237 map.push_back(iEdge_20[n]);
240 for(
unsigned int n=0; n<iEdge_03.size(); ++n){
242 map.push_back(iEdge_03[n]);
245 for(
unsigned int n=0; n<iEdge_13.size(); ++n){
247 map.push_back(iEdge_13[n]);
250 for(
unsigned int n=0; n<iEdge_23.size(); ++n){
252 map.push_back(iEdge_23[n]);
255 for(
unsigned int n=0; n<iFace_012.size(); ++n){
257 map.push_back(iFace_012[n]);
260 for(
unsigned int n=0; n<iFace_013.size(); ++n){
262 map.push_back(iFace_013[n]);
265 for(
unsigned int n=0; n<iFace_123.size(); ++n){
267 map.push_back(iFace_123[n]);
270 for(
unsigned int n=0; n<iFace_203.size(); ++n){
272 map.push_back(iFace_203[n]);
275 for(
unsigned int n=0; n<interiorVolumePoints.size(); ++n){
277 map.push_back(interiorVolumePoints[n]);
285 for(
unsigned int index=0; index<map.size(); ++index){
287 points[0][index] =
m_points[0][index];
288 points[1][index] =
m_points[1][index];
289 points[2][index] =
m_points[2][index];
293 for(
unsigned int index=0; index<map.size(); ++index){
295 m_points[0][index] = points[0][map[index]];
296 m_points[1][index] = points[1][map[index]];
297 m_points[2][index] = points[2][map[index]];
336 for(
int i = 0; i < rows; ++i ) {
337 for(
int j = 0; j < cols; ++j ) {
338 interp[j + i*cols] = interMat(i,j);
371 returnval->Initialize();
virtual void CalculatePoints()
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void CalculateInterpMatrix(const Array< OneD, const NekDouble > &xi, const Array< OneD, const NekDouble > &yi, const Array< OneD, const NekDouble > &zi, Array< OneD, NekDouble > &interp)
NekMatrix< NekDouble > GetTetInterpolationMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
MatrixSharedPtrType m_derivmatrix[3]
void CalculateDerivMatrix()
Array< OneD, DataType > m_points[3]
NekVector< NekDouble > MakeTetWeights(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z)
static boost::shared_ptr< PointsBaseType > Create(const PointsKey &key)
void NodalPointReorder3d()
unsigned int GetNumPoints() const
Array< OneD, DataType > m_weights
Defines a specification for a set of points.
Points< NekDouble >::MatrixSharedPtrType GetTetZDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
virtual void CalculateDerivMatrix()
Points< NekDouble >::MatrixSharedPtrType GetTetXDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
unsigned int GetRows() const
unsigned int GetTotNumPoints() const
Array< OneD, DataType > & GetPtr()
Points< NekDouble >::MatrixSharedPtrType GetTetYDerivativeMatrix(const NekVector< NekDouble > &x, const NekVector< NekDouble > &y, const NekVector< NekDouble > &z, const NekVector< NekDouble > &xi, const NekVector< NekDouble > &yi, const NekVector< NekDouble > &zi)
virtual void CalculateWeights()