Commit d5b1ccab authored by Michael Turner's avatar Michael Turner
Browse files

prisms not working correctly at higher orders

parent bbed186c
......@@ -37,7 +37,6 @@
#include <LibUtilities/Foundations/Points.h>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicConst/NektarUnivConsts.hpp>
#include <LibUtilities/Foundations/NodalTriSPIData.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Foundations/NodalUtil.h>
......@@ -51,28 +50,27 @@ namespace Nektar
unsigned int numPoints = GetNumPoints();
PointsKey e(numPoints,eGaussLobattoLegendre);
Array<OneD, NekDouble> gll;
PointsManager()[e]->GetPoints(gll);
PointsManager()[e]->GetPoints(e0);
ew = PointsManager()[e]->GetW();
PointsKey t(numPoints,eNodalTriSPI);
PointsManager()[t]->GetPoints(t0,t1);
tw = PointsManager()[t]->GetW();
numtri = tw.num_elements();
for(int i = 0; i < 3; i++)
{
m_points[i] = Array<OneD, DataType>(NodalTriSPINPTS[numPoints-2]*numPoints);
m_points[i] = Array<OneD, DataType>(numtri*numPoints);
}
for(int j = 0, ct = 0; j < numPoints; j++)
{
int index=0;
for(unsigned int i=0; i < numPoints-2; ++i)
{
index += NodalTriSPINPTS[i];
}
for(int i = 0; i < NodalTriSPINPTS[numPoints-2]; i++, ct++)
for(int i = 0; i < numtri; i++, ct++)
{
//need to flip y and z because of quad orientation
m_points[0][ct] = NodalTriSPIData[index][0];
m_points[1][ct] = gll[j];
m_points[2][ct] = NodalTriSPIData[index][1];
index++;
m_points[0][ct] = t0[i];
m_points[1][ct] = e0[j];
m_points[2][ct] = t1[i];
}
}
// exit(0);
......@@ -82,22 +80,13 @@ namespace Nektar
{
unsigned int numPoints = GetNumPoints();
m_weights = Array<OneD, DataType>(NodalTriSPINPTS[numPoints-2]*numPoints);
PointsKey e(numPoints,eGaussLobattoLegendre);
Array<OneD, NekDouble> gll = PointsManager()[e]->GetW();
m_weights = Array<OneD, DataType>(numtri*numPoints);
for(int j = 0, ct = 0; j < numPoints; j++)
{
int index=0;
for(unsigned int i=0; i < numPoints-2; ++i)
{
index += NodalTriSPINPTS[i];
}
for(int i = 0; i < NodalTriSPINPTS[numPoints-2]; i++, ct++)
for(int i = 0; i < numtri; i++, ct++)
{
m_weights[ct] = NodalTriSPIData[index][2] * gll[j];
index++;
m_weights[ct] = tw[i] * ew[j];
}
}
}
......
......@@ -68,6 +68,9 @@ namespace Nektar
{
}
Array<OneD, NekDouble> t0,t1,tw,e0,ew;
int numtri;
void CalculatePoints();
void CalculateWeights();
void CalculateDerivMatrix();
......
......@@ -153,6 +153,10 @@ void Edge::MakeOrder(int order,
loc[2] = m_edgeNodes[i]->m_z;
NekDouble t = CADCurve->loct(loc);
m_edgeNodes[i]->SetCADCurve(CADCurveId,CADCurve,t);
loc = CADCurve->P(t);
m_edgeNodes[i]->m_x = loc[0];
m_edgeNodes[i]->m_y = loc[1];
m_edgeNodes[i]->m_z = loc[2];
std::vector<CADSurfSharedPtr> s = CADCurve->GetAdjSurf();
for(int j = 0; j < s.size(); j++)
......@@ -171,6 +175,10 @@ void Edge::MakeOrder(int order,
loc[1] = m_edgeNodes[i]->m_y;
loc[2] = m_edgeNodes[i]->m_z;
Array<OneD, NekDouble> uv = CADSurf->locuv(loc);
loc = CADSurf->P(uv);
m_edgeNodes[i]->m_x = loc[0];
m_edgeNodes[i]->m_y = loc[1];
m_edgeNodes[i]->m_z = loc[2];
m_edgeNodes[i]->SetCADSurf(CADSurfId,CADSurf,uv);
}
}
......
......@@ -260,6 +260,10 @@ void Face::MakeOrder(int order,
loc[1] = m_faceNodes[i]->m_y;
loc[2] = m_faceNodes[i]->m_z;
Array<OneD, NekDouble> uv = CADSurf->locuv(loc);
loc = CADSurf->P(uv);
m_faceNodes[i]->m_x = loc[0];
m_faceNodes[i]->m_y = loc[1];
m_faceNodes[i]->m_z = loc[2];
m_faceNodes[i]->SetCADSurf(CADSurfId,CADSurf,uv);
}
}
......
......@@ -217,7 +217,7 @@ void ProcessVarOpti::BuildDerivUtil()
derivUtil[st]->ptsLow = m_mesh->m_nummode*m_mesh->m_nummode*(m_mesh->m_nummode+1)/2;
LibUtilities::PointsKey pkey1(m_mesh->m_nummode,
LibUtilities::eNodalPrismElec);
LibUtilities::PointsKey pkey2(m_mesh->m_nummode+5,
LibUtilities::PointsKey pkey2(m_mesh->m_nummode+4,
LibUtilities::eNodalPrismSPI);
Array<OneD, NekDouble> u1, v1, u2, v2, w1, w2;
LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1, w1);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment