Commit 4a3d0952 authored by Dave Moxey's avatar Dave Moxey

Fix some silly errors in the derivative modes setup, move things so that CAD is not required

parent f3d6a5de
......@@ -470,8 +470,8 @@ public:
NEKMESHUTILS_EXPORT virtual void GetCurvedNodes(
std::vector<NodeSharedPtr> &nodeList) const
{
cerr << "WARNING: Unsupported curvature for a " << m_vertex.size()
<< "-vertex element is not yet implemented." << endl;
std::cerr << "WARNING: Unsupported curvature for a " << m_vertex.size()
<< "-vertex element is not yet implemented." << std::endl;
// Node orderings are different for different elements.
// Triangle
if (m_vertex.size() == 2)
......
......@@ -354,7 +354,21 @@ public:
}
return ret;
}
#else
int GetNumCadCurve()
{
return 0;
}
int GetNumCADSurf()
{
return 0;
}
std::vector<std::string> GetCADString()
{
return std::vector<std::string>();
}
#endif
/// ID of node.
......
......@@ -23,6 +23,10 @@ SET(NekMeshHeaders
ProcessModules/ProcessScalar.h
ProcessModules/ProcessSpherigon.h
ProcessModules/ProcessTetSplit.h
ProcessModules/ProcessOptiExtract.h
ProcessModules/ProcessVarOpti/ProcessVarOpti.h
ProcessModules/ProcessVarOpti/NodeOpti.h
ProcessModules/ProcessVarOpti/ElUtil.h
)
SET(NekMeshSources
......@@ -52,6 +56,10 @@ SET(NekMeshSources
ProcessModules/ProcessSpherigon.cpp
ProcessModules/ProcessTetSplit.cpp
ProcessModules/ProcessOptiExtract.cpp
ProcessModules/ProcessVarOpti/ProcessVarOpti.cpp
ProcessModules/ProcessVarOpti/PreProcessing.cpp
ProcessModules/ProcessVarOpti/NodeOpti.cpp
ProcessModules/ProcessVarOpti/ElUtil.cpp
)
IF (NEKTAR_USE_CCM)
......@@ -66,17 +74,11 @@ ENDIF (NEKTAR_USE_VTK)
IF (NEKTAR_USE_MESHGEN)
SET(NekMeshHeaders ${NekMeshHeaders}
InputModules/InputCAD.h
ProcessModules/ProcessOptiExtract.h
ProcessModules/ProcessVarOpti/ProcessVarOpti.h
ProcessModules/ProcessVarOpti/NodeOpti.h
ProcessModules/ProcessVarOpti/ElUtil.h)
ProcessModules/ProcessVarOpti/NodeOptiCAD.h
InputModules/InputCAD.h)
SET(NekMeshSources ${NekMeshSources}
InputModules/InputCAD.cpp
ProcessModules/ProcessVarOpti/ProcessVarOpti.cpp
ProcessModules/ProcessVarOpti/PreProcessing.cpp
ProcessModules/ProcessVarOpti/NodeOpti.cpp
ProcessModules/ProcessVarOpti/ElUtil.cpp)
ProcessModules/ProcessVarOpti/NodeOptiCAD.cpp
InputModules/InputCAD.cpp)
ENDIF (NEKTAR_USE_MESHGEN)
IF (NEKTAR_USE_ANN)
......
......@@ -1028,7 +1028,7 @@ void OutputNekpp::WriteXmlCAD(TiXmlElement *pRoot)
}
}
}
{
EdgeSet::iterator it;
for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end();
......@@ -1100,7 +1100,10 @@ void OutputNekpp::WriteXmlCAD(TiXmlElement *pRoot)
}
}
pRoot->LinkEndChild(cad);
if (cad->FirstChild())
{
pRoot->LinkEndChild(cad);
}
}
void OutputNekpp::WriteXmlCADId(TiXmlElement *pRoot)
......
......@@ -139,58 +139,36 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
}
else if(m_el->GetConf().m_e == LibUtilities::eTriangle)
{
/*LibUtilities::PointsKey pkey(m_mode,
LibUtilities::eNodalTriElec);
Array<OneD, NekDouble> u, v;
LibUtilities::PointsManager()[pkey]->GetPoints(u, v);
Array<OneD, NekDouble> xc(chi->GetTotPoints());
Array<OneD, NekDouble> yc(chi->GetTotPoints());
Array<OneD, NekDouble> coeffs0 = geom->GetCoeffs(0);
Array<OneD, NekDouble> coeffs1 = geom->GetCoeffs(1);
chi->BwdTrans(coeffs0,xc);
chi->BwdTrans(coeffs1,yc);
NekVector<NekDouble> X(ptsHelp->ptsLow),Y(ptsHelp->ptsLow);
for(int j = 0; j < u.num_elements(); j++)
{
Array<OneD, NekDouble> xp(2);
xp[0] = u[j];
xp[1] = v[j];
DNekMat J(2,2,0.0);
J(0,0) = (*nodes[1][0] - *nodes[0][0]);
J(1,0) = (*nodes[1][1] - *nodes[0][1]);
J(0,1) = (*nodes[2][0] - *nodes[0][0]);
J(1,1) = (*nodes[2][1] - *nodes[0][1]);
X(j) = chi->PhysEvaluate(xp, xc);
Y(j) = chi->PhysEvaluate(xp, yc);
}
J.Invert();
NekVector<NekDouble> x1i(ptsHelp->ptsHigh),y1i(ptsHelp->ptsHigh),
x2i(ptsHelp->ptsHigh),y2i(ptsHelp->ptsHigh);
DNekMat R(2,2,0.0);
R(0,0) = 2.0;
R(1,1) = 2.0;
x1i = derivUtil->VdmD[0]*X;
y1i = derivUtil->VdmD[0]*Y;
x2i = derivUtil->VdmD[1]*X;
y2i = derivUtil->VdmD[1]*Y;
J = J * R;
for(int i = 0 ; i < ptsHelp->ptsHigh; i++)
for(int i = 0 ; i < derivUtil->ptsHigh; i++)
{
DNekMat dxdz(2,2,1.0,eFULL);
dxdz(0,0) = x1i(i);
dxdz(0,1) = x2i(i);
dxdz(1,0) = y1i(i);
dxdz(1,1) = y2i(i);
Array<OneD, NekDouble> r(10,0.0);
r[9] = dxdz(0,0)*dxdz(1,1)-dxdz(1,0)*dxdz(0,1);
dxdz.Invert();
Array<OneD, NekDouble> r(10,0.0); //store det in 10th entry
r[0] = dxdz(0,0);
r[1] = dxdz(1,0);
r[3] = dxdz(0,1);
r[4] = dxdz(1,1);
r[9] = 1.0 / (J(0,0) * J(1,1) - J(0,1) * J(1,0));
r[0] = J(0,0);
r[1] = J(1,0);
r[2] = 0.0;
r[3] = J(0,1);
r[4] = J(1,1);
r[5] = 0.0;
r[6] = 0.0;
r[7] = 0.0;
r[8] = 0.0;
ret.push_back(r);
}*/
}
}
else if(m_el->GetConf().m_e == LibUtilities::eTetrahedron)
{
......@@ -396,7 +374,7 @@ void ElUtil::Evaluate()
}
//delta = minEdge / m_el->GetConf().m_order / 500.0;
delta = 1e-4;
delta = 1e-6;
minJac = mn;
}
......
......@@ -43,7 +43,12 @@ namespace Nektar
namespace Utilities
{
boost::mutex mtx;
NodeOptiFactory &GetNodeOptiFactory()
{
typedef Loki::SingletonHolder<NodeOptiFactory, Loki::CreateUsingNew,
Loki::NoDestroy> Type;
return Type::Instance();
}
void NodeOpti::CalcDX()
{
......@@ -65,6 +70,9 @@ void NodeOpti::CalcMinJac()
}
}
int NodeOpti2D2D::m_type = GetNodeOptiFactory().RegisterCreatorFunction(
22, NodeOpti2D2D::create, "2D2D");
void NodeOpti2D2D::Optimise()
{
CalcDX();
......@@ -110,6 +118,42 @@ void NodeOpti2D2D::Optimise()
}
}
Array<OneD, NekDouble> NodeOpti2D2D::GetGrad()
{
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
vector<NekDouble> w(9);
for(int i = 0; i < 7; i++)
{
node->m_x = xc + dir[i][0] * dx;
node->m_y = yc + dir[i][1] * dx;
w[i] = GetFunctional<2>();
}
node->m_x = xc;
node->m_y = yc;
Array<OneD, NekDouble> ret(5,0.0);
//ret[0] d/dx
//ret[1] d/dy
//ret[3] d2/dx2
//ret[4] d2/dy2
//ret[5] d2/dxdy
ret[0] = (w[1] - w[4]) / 2.0 / dx;
ret[1] = (w[3] - w[6]) / 2.0 / dx;
ret[2] = (w[1] + w[4] - 2.0*w[0]) / dx / dx;
ret[3] = (w[3] + w[6] - 2.0*w[0]) / dx / dx;
ret[4] = (w[2] - w[1] - w[3] + 2.0*w[0] - w[4] - w[6] + w[5]) / 2.0 / dx / dx;
return ret;
}
int NodeOpti3D3D::m_type = GetNodeOptiFactory().RegisterCreatorFunction(
33, NodeOpti3D3D::create, "3D3D");
void NodeOpti3D3D::Optimise()
{
......@@ -136,14 +180,20 @@ void NodeOpti3D3D::Optimise()
-G[6]*(G[6]*G[5]-G[7]*G[8])
+G[7]*(G[6]*G[8]-G[7]*G[4]);
delX = G[0]*(G[4]*G[5]-G[8]*G[8]) + G[1]*(G[7]*G[8]-G[6]*G[5]) + G[2]*(G[6]*G[8]-G[7]*G[4]);
delX = G[0]*(G[4]*G[5]-G[8]*G[8]) +
G[1]*(G[7]*G[8]-G[6]*G[5]) +
G[2]*(G[6]*G[8]-G[7]*G[4]);
delY = G[0]*(G[8]*G[7]-G[6]*G[5]) +
G[1]*(G[3]*G[5]-G[7]*G[7]) +
G[2]*(G[6]*G[7]-G[3]*G[8]);
delZ = G[0]*(G[6]*G[8]-G[4]*G[7]) +
G[1]*(G[6]*G[7]-G[3]*G[8]) +
G[2]*(G[3]*G[4]-G[6]*G[6]);
delX /= det;
delY = G[0]*(G[8]*G[7]-G[6]*G[5]) + G[1]*(G[3]*G[5]-G[7]*G[7]) + G[2]*(G[6]*G[7]-G[3]*G[8]);
delY /= det;
delZ = G[0]*(G[6]*G[8]-G[4]*G[7]) + G[1]*(G[6]*G[7]-G[3]*G[8]) + G[2]*(G[3]*G[4]-G[6]*G[6]);
delZ /= det;
bool found = false;
while(alpha > 1e-10)
{
......@@ -166,7 +216,7 @@ void NodeOpti3D3D::Optimise()
node->m_x = xc;
node->m_y = yc;
node->m_z = zc;
functional = currentW;
//functional = currentW;
//cout << "warning: had to reset node" << endl;
}
mtx.lock();
......@@ -176,288 +226,6 @@ void NodeOpti3D3D::Optimise()
}
}
void NodeOpti1D3D::Optimise()
{
CalcDX();
CalcMinJac();
Array<OneD, NekDouble> G = GetGrad();
if(sqrt(G[0]*G[0]) > 1e-10)
{
//needs to optimise
NekDouble tc = node->GetCADCurveInfo(curve->GetId());
NekDouble currentW = GetFunctional<3>();
NekDouble functional;
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
NekDouble zc = node->m_z;
NekDouble alpha = 1.0;
NekDouble delT;
NekDouble nt;
Array<OneD, NekDouble> p;
delT = G[0] / G[1];
Array<OneD, NekDouble> bd = curve->Bounds();
bool found = false;
while(alpha > 1e-10)
{
nt = tc - alpha * delT;
if(nt < bd[0] || nt > bd[1])
{
alpha /= 2.0;
continue;
}
p = curve->P(nt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
functional = GetFunctional<3>();
if(functional < currentW)
{
found = true;
break;
}
alpha /= 2.0;
}
if(!found)
{
//reset the node
nt = tc;
p = curve->P(nt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
functional = currentW;
// cout << "warning: had to reset node" << endl;
}
else
{
node->MoveCurve(p,curve->GetId(),nt);
}
mtx.lock();
res->val = max(sqrt((node->m_x-xc)*(node->m_x-xc)+(node->m_y-yc)*(node->m_y-yc)+
(node->m_z-zc)*(node->m_z-zc)),res->val);
mtx.unlock();
}
}
void NodeOpti2D3D::Optimise()
{
CalcDX();
CalcMinJac();
Array<OneD, NekDouble> G = GetGrad();
if(sqrt(G[0]*G[0] + G[1]*G[1]) > 1e-10)
{
//needs to optimise
Array<OneD, NekDouble> uvc = node->GetCADSurfInfo(surf->GetId());
NekDouble currentW = GetFunctional<3>();
NekDouble functional;
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
NekDouble zc = node->m_z;
NekDouble alpha = 1.0;
Array<OneD, NekDouble> uvt(2);
Array<OneD, NekDouble> p;
NekDouble delU = 1.0/(G[2]*G[3]-G[4]*G[4])*(G[3]*G[0] - G[4]*G[1]);
NekDouble delV = 1.0/(G[2]*G[3]-G[4]*G[4])*(G[2]*G[1] - G[4]*G[0]);
Array<OneD, NekDouble> bd = surf->GetBounds();
bool found = false;
while(alpha > 1e-10)
{
uvt[0] = uvc[0] - alpha * delU;
uvt[1] = uvc[1] - alpha * delV;
if(uvt[0] < bd[0] || uvt[0] > bd[1] ||
uvt[1] < bd[2] || uvt[1] > bd[3])
{
alpha /= 2.0;
continue;
}
p = surf->P(uvt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
functional = GetFunctional<3>();
if(functional < currentW)
{
found = true;
break;
}
alpha /= 2.0;
}
if(!found)
{
//reset the node
p = surf->P(uvc);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
functional = currentW;
// cout << "warning: had to reset node" << endl;
}
else
{
node->Move(p,surf->GetId(),uvt);
}
mtx.lock();
res->val = max(sqrt((node->m_x-xc)*(node->m_x-xc)+(node->m_y-yc)*(node->m_y-yc)+
(node->m_z-zc)*(node->m_z-zc)),res->val);
mtx.unlock();
}
}
NekDouble dir[13][3] = {{ 0.0, 0.0, 0.0 }, // 0 (x , y , z )
{ 1.0, 0.0, 0.0 }, // 1 (x+dx, y , z )
{ 1.0, 1.0, 0.0 }, // 2 (x+dx, y+dy, z )
{ 0.0, 1.0, 0.0 }, // 3 (x , y+dy, z )
{ -1.0, 0.0, 0.0 }, // 4 (x-dx, y , z )
{ -1.0, -1.0, 0.0 }, // 5 (x-dx, y-dy, z )
{ 0.0, -1.0, 0.0 }, // 6 (x , y-dy, z )
{ -1.0, 0.0, -1.0 }, // 7 (x-dx, y , z-dz)
{ 0.0, 0.0, -1.0 }, // 8 (x , y , z-dz)
{ 0.0, 0.0, 1.0 }, // 9 (x , y , z+dz)
{ 1.0, 0.0, 1.0 }, // 10 (x+dx, y , z+dz)
{ 0.0, 1.0, 1.0 }, // 11 (x , y+dy, z+dz)
{ 0.0, -1.0, -1.0 }}; // 12 (x , y-dy, z-dz)
Array<OneD, NekDouble> NodeOpti1D3D::GetGrad()
{
NekDouble tc = node->GetCADCurveInfo(curve->GetId());
Array<OneD, NekDouble> d1 = curve->D1(tc);
NekDouble dr = sqrt(d1[3]*d1[3] + d1[4]*d1[4] + d1[5]*d1[5]);
NekDouble dt = dx / dr;
vector<NekDouble> w(3);
w[0] = GetFunctional<3>();
NekDouble nt = tc + dt;
Array<OneD, NekDouble> p = curve->P(nt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
w[1] = GetFunctional<3>();
nt = tc - dt;
p = curve->P(nt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
w[2] = GetFunctional<3>();
nt = tc;
p = curve->P(nt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
Array<OneD, NekDouble> ret(2,0.0);
ret[0] = (w[1] - w[2]) / 2.0 / dt;
ret[1] = (w[1] + w[2] - 2.0*w[0]) / dt / dt;
return ret;
}
Array<OneD, NekDouble> NodeOpti2D3D::GetGrad()
{
Array<OneD, NekDouble> uvc = node->GetCADSurfInfo(surf->GetId());
Array<OneD, NekDouble> d1 = surf->D1(uvc);
NekDouble dru = sqrt(d1[3]*d1[3] + d1[4]*d1[4] + d1[5]*d1[5]);
NekDouble drv = sqrt(d1[6]*d1[6] + d1[7]*d1[7] + d1[8]*d1[8]);
NekDouble du = dx / dru;
NekDouble dv = dx / drv;
vector<NekDouble> w(7);
for(int i = 0; i < 7; i++)
{
Array<OneD, NekDouble> uvt(2);
uvt[0] = uvc[0] + dir[i][0] * du;
uvt[1] = uvc[1] + dir[i][1] * dv;
Array<OneD, NekDouble> p = surf->P(uvt);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
w[i] = GetFunctional<3>();
}
Array<OneD, NekDouble> p = surf->P(uvc);
node->m_x = p[0];
node->m_y = p[1];
node->m_z = p[2];
Array<OneD, NekDouble> ret(5,0.0);
//ret[0] d/dx
//ret[1] d/dy
//ret[2] d2/dx2
//ret[3] d2/dy2
//ret[4] d2/dxdy
ret[0] = (w[1] - w[4]) / 2.0 / du;
ret[1] = (w[3] - w[6]) / 2.0 / dv;
ret[2] = (w[1] + w[4] - 2.0*w[0]) / du / du;
ret[3] = (w[3] + w[6] - 2.0*w[0]) / dv / dv;
ret[4] = (w[2] - w[1] - w[3] + 2.0*w[0] - w[4] - w[6] + w[5]) / 2.0 / du / dv;
return ret;
}
Array<OneD, NekDouble> NodeOpti2D2D::GetGrad()
{
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
vector<NekDouble> w(9);
for(int i = 0; i < 7; i++)
{
node->m_x = xc + dir[i][0] * dx;
node->m_y = yc + dir[i][1] * dx;
w[i] = GetFunctional<2>();
}
node->m_x = xc;
node->m_y = yc;
Array<OneD, NekDouble> ret(5,0.0);
//ret[0] d/dx
//ret[1] d/dy
//ret[3] d2/dx2
//ret[4] d2/dy2
//ret[5] d2/dxdy
ret[0] = (w[1] - w[4]) / 2.0 / dx;
ret[1] = (w[3] - w[6]) / 2.0 / dx;
ret[2] = (w[1] + w[4] - 2.0*w[0]) / dx / dx;
ret[3] = (w[3] + w[6] - 2.0*w[0]) / dx / dx;
ret[4] = (w[2] - w[1] - w[3] + 2.0*w[0] - w[4] - w[6] + w[5]) / 2.0 / dx / dx;
return ret;
}
Array<OneD, NekDouble> NodeOpti3D3D::