Commit 39a9d48a authored by Michael Turner's avatar Michael Turner
Browse files

prisms working for the linear case

parent b369ab5b
......@@ -138,7 +138,9 @@ SET(FoundationHeaders
./Foundations/NodalTriSPIData.h
./Foundations/NodalTriSPI.h
./Foundations/NodalTetSPI.h
./Foundations/NodalPrismSPI.h
./Foundations/NodalTetSPIData.h
./Foundations/NodalPrismSPIData.h
./Foundations/NodalQuadEvenlySpaced.h
./Foundations/NodalQuadElec.h
./Foundations/NodalUtil.h
......@@ -167,6 +169,7 @@ SET(FoundationSources
./Foundations/NodalTriFekete.cpp
./Foundations/NodalTriSPI.cpp
./Foundations/NodalTetSPI.cpp
./Foundations/NodalPrismSPI.cpp
./Foundations/NodalQuadEvenlySpaced.cpp
./Foundations/NodalQuadElec.cpp
./Foundations/NodalUtil.cpp
......
......@@ -100,7 +100,8 @@ namespace Nektar
"NodalTetElec",
"NodalTetSPI",
"NodalPrismEvenlySpaced"
"NodalPrismElec"
"NodalPrismElec",
"NodalPrismSPI"
};
} // end of namespace
......
......@@ -46,6 +46,7 @@
#include <LibUtilities/Foundations/NodalTriEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalTetEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalTetSPI.h>
#include <LibUtilities/Foundations/NodalPrismSPI.h>
#include <LibUtilities/Foundations/NodalPrismEvenlySpaced.h>
#include <LibUtilities/Foundations/NodalPrismElec.h>
#include <LibUtilities/Foundations/NodalQuadEvenlySpaced.h>
......@@ -95,6 +96,7 @@ namespace Nektar
const bool NodalTetEveInited = PointsManager().RegisterCreator(PointsKey(0, eNodalTetEvenlySpaced), NodalTetEvenlySpaced::Create);
const bool NodalPrismEveInited = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismEvenlySpaced), NodalPrismEvenlySpaced::Create);
const bool NodalPrismElecInited = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismElec), NodalPrismElec::Create);
const bool NodalPrismSPIInited = PointsManager().RegisterCreator(PointsKey(0, eNodalPrismSPI), NodalPrismSPI::Create);
const bool Basis_Inited = BasisManager().RegisterGlobalCreator(Basis::Create);
};
......
......@@ -249,9 +249,7 @@ namespace Nektar
for(unsigned int index=0; index<map.size(); ++index){
m_points[0][index] = points[0][map[index]];
m_points[1][index] = points[1][map[index]];
cout << m_points[0][index] << " " << m_points[1][index] << endl;
}
exit(-1);
}
......
......@@ -164,6 +164,7 @@ namespace Nektar
case eNodalTetEvenlySpaced:
case eNodalPrismEvenlySpaced:
case eNodalPrismElec:
case eNodalPrismSPI:
dimpoints = 3;
break;
......@@ -200,6 +201,9 @@ namespace Nektar
case eNodalPrismElec:
totpoints = m_numpoints*m_numpoints*(m_numpoints+1)/2;
break;
case eNodalPrismSPI: ASSERTL0(false,"this method cannot be implemented");
break;
break;
case eNodalQuadEvenlySpaced:
case eNodalQuadElec:
......
......@@ -78,6 +78,7 @@ namespace Nektar
eNodalTetSPI, //!< 3D Nodal Symmetric positive internal tet (Whitherden, Vincent)
eNodalPrismEvenlySpaced, //!< 3D Evenly-spaced points on a Prism
eNodalPrismElec, //!< 3D electrostatically spaced points on a Prism
eNodalPrismSPI, //!< 3D prism SPI
SIZE_PointsType //!< Length of enum list
};
......
......@@ -479,19 +479,6 @@ void Prism::GetCurvedNodes(std::vector<NodeSharedPtr> &nodeList) const
fcid.push_back(m_face[i]->m_vertexList[2]->m_id);
fcid.push_back(m_face[i]->m_vertexList[3]->m_id);
for(int j = 0; j < fcid.size(); j++)
{
cout << m_face[i]->m_vertexList[j]->m_x << " " << m_face[i]->m_vertexList[j]->m_y << " " << m_face[i]->m_vertexList[j]->m_z << endl;
}
for(int j = 0; j < m_face[i]->m_faceNodes.size(); j++)
{
cout << m_face[i]->m_faceNodes[j]->m_x << " " << m_face[i]->m_faceNodes[j]->m_y << " " << m_face[i]->m_faceNodes[j]->m_z << endl;
}
exit(-1);
HOQuadrilateral<NodeSharedPtr> hoq(fcid, m_face[i]->m_faceNodes);
hoq.Align(ts[i]);
......@@ -506,12 +493,6 @@ void Prism::GetCurvedNodes(std::vector<NodeSharedPtr> &nodeList) const
std::copy(m_volumeNodes.begin(),
m_volumeNodes.end(),
nodeList.begin() + k);
for(int i = 0; i < nodeList.size(); i++)
{
cout << nodeList[i]->m_x << " " << nodeList[i]->m_y << " " << nodeList[i]->m_z << endl;
}
exit(-1);
}
/**
......
......@@ -78,25 +78,10 @@ ElUtil::ElUtil(ElementSharedPtr e, DerivUtilSharedPtr d,
vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
{
//need to make ideal element out of old element
ElmtConfig ec = m_el->GetConf();
ec.m_order = 1;
ec.m_faceNodes = false;
ec.m_volumeNodes = false;
ec.m_reorient = false;
ElementSharedPtr E = GetElementFactory().CreateInstance(
ec.m_e, ec, m_el->GetVertexList(),
m_el->GetTagList());
mtx2.lock();
SpatialDomains::GeometrySharedPtr geom = E->GetGeom(m_dim);
geom->FillGeom();
StdRegions::StdExpansionSharedPtr chi = geom->GetXmap();
mtx2.unlock();
vector<Array<OneD, NekDouble> > ret;
LibUtilities::ShapeType st = m_el->GetConf().m_e;
if(m_el->GetConf().m_e == LibUtilities::eQuadrilateral)
{
ASSERTL0(false,"Not coded");
......@@ -216,60 +201,76 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
}
else if(m_el->GetConf().m_e == LibUtilities::ePrism)
{
vector<Array<OneD, NekDouble> > xyz;
for(int i = 0; i < geom->GetNumVerts(); i++)
// r,s,t are the coordinates of the SPI points
Array<OneD, NekDouble> r = derivUtil->ptx;
Array<OneD, NekDouble> s = derivUtil->pty;
Array<OneD, NekDouble> t = derivUtil->ptz;
vector<Array<OneD, NekDouble> > xyz(6);
vector<NodeSharedPtr> ns = m_el->GetVertexList();
for(int i = 0; i < 6; i++)
{
Array<OneD, NekDouble> loc(3);
SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
p->GetCoords(loc);
xyz.push_back(loc);
Array<OneD, NekDouble> x(3);
x[0] = ns[i]->m_x;
x[1] = ns[i]->m_y;
x[2] = ns[i]->m_z;
xyz[i] = x;
}
Array<OneD, const LibUtilities::BasisSharedPtr> b = chi->GetBase();
Array<OneD, NekDouble> eta1 = b[0]->GetZ();
Array<OneD, NekDouble> eta2 = b[1]->GetZ();
Array<OneD, NekDouble> eta3 = b[2]->GetZ();
for(int k = 0; k < b[2]->GetNumPoints(); k++)
for (int i = 0; i < r.num_elements(); ++i)
{
for(int j = 0; j < b[1]->GetNumPoints(); j++)
{
for(int i = 0; i < b[0]->GetNumPoints(); i++)
{
NekDouble xi1 = 0.5*(1+eta1[i])*(1-eta3[k])-1.0;
NekDouble a1 = 0.5*(1-xi1), a2 = 0.5*(1+xi1);
NekDouble b1 = 0.5*(1-eta2[j]), b2 = 0.5*(1+eta2[j]);
NekDouble c1 = 0.5*(1-eta3[k]), c2 = 0.5*(1+eta3[k]);
DNekMat dxdz(3,3,1.0,eFULL);
dxdz(0,0) = 0.5*(-b1*xyz[0][0] + b1*xyz[1][0] + b2*xyz[2][0] - b2*xyz[3][0]);
dxdz(1,0) = 0.5*(-b1*xyz[0][1] + b1*xyz[1][1] + b2*xyz[2][1] - b2*xyz[3][1]);
dxdz(2,0) = 0.5*(-b1*xyz[0][2] + b1*xyz[1][2] + b2*xyz[2][2] - b2*xyz[3][2]);
dxdz(0,1) = 0.5*((a2-c1)*xyz[0][0] - a2*xyz[1][0] + a2*xyz[2][0] + (c1-a2)*xyz[3][0] - c2*xyz[4][0] + c2*xyz[5][0]);
dxdz(1,1) = 0.5*((a2-c1)*xyz[0][1] - a2*xyz[1][1] + a2*xyz[2][1] + (c1-a2)*xyz[3][1] - c2*xyz[4][1] + c2*xyz[5][1]);
dxdz(2,1) = 0.5*((a2-c1)*xyz[0][2] - a2*xyz[1][2] + a2*xyz[2][2] + (c1-a2)*xyz[3][2] - c2*xyz[4][2] + c2*xyz[5][2]);
dxdz(0,2) = 0.5*(-b1*xyz[0][0] - b2*xyz[3][0] + b1*xyz[4][0] + b2*xyz[5][0]);
dxdz(1,2) = 0.5*(-b1*xyz[0][1] - b2*xyz[3][1] + b1*xyz[4][1] + b2*xyz[5][1]);
dxdz(2,2) = 0.5*(-b1*xyz[0][2] - b2*xyz[3][2] + b1*xyz[4][2] + b2*xyz[5][2]);
dxdz.Invert();
Array<OneD, NekDouble> r(9,0.0);
r[0] = dxdz(0,0);
r[1] = dxdz(1,0);
r[3] = dxdz(0,1);
r[4] = dxdz(1,1);
r[2] = dxdz(2,0);
r[5] = dxdz(2,1);
r[6] = dxdz(0,2);
r[7] = dxdz(1,2);
r[8] = dxdz(2,2);
ret.push_back(r);
}
}
NekDouble a2 = 0.5 * (1 + r[i]);
NekDouble b1 = 0.5 * (1 - s[i]);
NekDouble b2 = 0.5 * (1 + s[i]);
NekDouble c1 = 0.5 * (1 - t[i]);
NekDouble c2 = 0.5 * (1 + t[i]);
DNekMat J(3, 3, 1.0, eFULL);
J(0, 0) = 0.5 * (-b1 * xyz[0][0] + b1 * xyz[1][0] +
b2 * xyz[2][0] - b2 * xyz[3][0]);
J(1, 0) = 0.5 * (-b1 * xyz[0][1] + b1 * xyz[1][1] +
b2 * xyz[2][1] - b2 * xyz[3][1]);
J(2, 0) = 0.5 * (-b1 * xyz[0][2] + b1 * xyz[1][2] +
b2 * xyz[2][2] - b2 * xyz[3][2]);
J(0, 1) = 0.5 * ((a2 - c1) * xyz[0][0] - a2 * xyz[1][0] +
a2 * xyz[2][0] + (c1 - a2) * xyz[3][0] -
c2 * xyz[4][0] + c2 * xyz[5][0]);
J(1, 1) = 0.5 * ((a2 - c1) * xyz[0][1] - a2 * xyz[1][1] +
a2 * xyz[2][1] + (c1 - a2) * xyz[3][1] -
c2 * xyz[4][1] + c2 * xyz[5][1]);
J(2, 1) = 0.5 * ((a2 - c1) * xyz[0][2] - a2 * xyz[1][2] +
a2 * xyz[2][2] + (c1 - a2) * xyz[3][2] -
c2 * xyz[4][2] + c2 * xyz[5][2]);
J(0, 2) = 0.5 * (-b1 * xyz[0][0] - b2 * xyz[3][0] +
b1 * xyz[4][0] + b2 * xyz[5][0]);
J(1, 2) = 0.5 * (-b1 * xyz[0][1] - b2 * xyz[3][1] +
b1 * xyz[4][1] + b2 * xyz[5][1]);
J(2, 2) = 0.5 * (-b1 * xyz[0][2] - b2 * xyz[3][2] +
b1 * xyz[4][2] + b2 * xyz[5][2]);
J.Invert();
Array<OneD, NekDouble> r(10,0.0); //store det in 10th entry
r[9] = 1.0/(J(0,0)*(J(1,1)*J(2,2)-J(2,1)*J(1,2))
-J(0,1)*(J(1,0)*J(2,2)-J(2,0)*J(1,2))
+J(0,2)*(J(1,0)*J(2,1)-J(2,0)*J(1,1)));
r[0] = J(0,0);
r[1] = J(1,0);
r[2] = J(2,0);
r[3] = J(0,1);
r[4] = J(1,1);
r[5] = J(2,1);
r[6] = J(0,2);
r[7] = J(1,2);
r[8] = J(2,2);
ret.push_back(r);
}
}
else
{
......
......@@ -150,7 +150,7 @@ void NodeOpti3D3D::Optimise()
NekDouble currentW = GetFunctional<3>();
NekDouble newVal;
if(G[0]*G[0] + G[1]*G[1] + G[2]*G[2] > gradTol())
{
//needs to optimise
......
......@@ -143,7 +143,7 @@ void ProcessVarOpti::BuildDerivUtil()
LibUtilities::PointsKey pkey1(m_mesh->m_nummode,
LibUtilities::eNodalPrismElec);
LibUtilities::PointsKey pkey2(m_mesh->m_nummode+4,
LibUtilities::eNodalPrismElec);
LibUtilities::eNodalPrismSPI);
Array<OneD, NekDouble> u1, v1, u2, v2, w1, w2;
LibUtilities::PointsManager()[pkey1]->GetPoints(u1, v1, w1);
LibUtilities::PointsManager()[pkey2]->GetPoints(u2, v2, w2);
......@@ -176,6 +176,10 @@ void ProcessVarOpti::BuildDerivUtil()
Array<OneD, NekDouble> qds = LibUtilities::PointsManager()[pkey2]->GetW();
NekVector<NekDouble> quadWi(qds);
derivUtil[st]->quadW = quadWi;
derivUtil[st]->ptx = u2;
derivUtil[st]->pty = v2;
derivUtil[st]->ptz = w2;
}
}
}
......
......@@ -286,14 +286,14 @@ void ProcessVarOpti::Process()
jobs[j] = optiNodes[i][j]->GetJob();
}
tm->SetNumWorkers(0);
tm->QueueJobs(jobs);
tm->SetNumWorkers(nThreads);
tm->Wait();
//for(int j = 0; j < jobs.size(); j++)
//{
// jobs[j]->Run();
//}
//tm->SetNumWorkers(0);
//tm->QueueJobs(jobs);
//tm->SetNumWorkers(nThreads);
//tm->Wait();
for(int j = 0; j < jobs.size(); j++)
{
jobs[j]->Run();
}
}
res->startInv = 0;
......
......@@ -55,6 +55,8 @@ struct DerivUtil
int ptsHigh;
int ptsLow;
Array<OneD, NekDouble> ptx, pty, ptz;
};
typedef boost::shared_ptr<DerivUtil> DerivUtilSharedPtr;
......
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