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

latest changes

parent 3f82d4d7
......@@ -55,14 +55,8 @@ public:
std::vector<NodeSharedPtr> pNodeList,
std::vector<int> pTagList)
{
ElementSharedPtr e = boost::shared_ptr<Element>(
return boost::shared_ptr<Element>(
new Quadrilateral(pConf, pNodeList, pTagList));
std::vector<EdgeSharedPtr> m_edges = e->GetEdgeList();
for (int i = 0; i < m_edges.size(); ++i)
{
m_edges[i]->m_elLink.push_back(std::pair<ElementSharedPtr, int>(e, i));
}
return e;
}
/// Element type
static LibUtilities::ShapeType m_type;
......
......@@ -83,6 +83,11 @@ public:
return m_idmap[in];
}
void CalcMap()
{
maps = MappingIdealToRef();
}
private:
std::vector<Array<OneD, NekDouble> > MappingIdealToRef();
......
......@@ -81,6 +81,18 @@ void NodeOpti::CalcMinJac()
}
}
void NodeOpti::ReCalcMaps()
{
std::map<LibUtilities::ShapeType,std::vector<ElUtilSharedPtr> >::iterator typeIt;
for(typeIt = data.begin(); typeIt != data.end(); typeIt++)
{
for(int i = 0; i < typeIt->second.size(); i++)
{
typeIt->second[i]->CalcMap();
}
}
}
int NodeOpti2D2D::m_type = GetNodeOptiFactory().RegisterCreatorFunction(
22, NodeOpti2D2D::create, "2D2D");
......@@ -98,11 +110,16 @@ void NodeOpti2D2D::Optimise()
NekDouble currentW = GetFunctional<2>(false,false);
NekDouble newVal;
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
vector<NekDouble> d;
for(int i = 0; i < 6; i++)
{
NekDouble xc = node->m_x + 1e-4*dir[i][0];
NekDouble yc = node->m_y + 1e-4*dir[i][1];
node->m_x = xc + 1e-4*dir[i][0];
node->m_y = yc + 1e-4*dir[i][1];
ReCalcMaps();
d.push_back(GetFunctional<2>(false,false));
}
......@@ -110,16 +127,14 @@ void NodeOpti2D2D::Optimise()
G = Array<OneD, NekDouble>(5);
G[0] = (d[1] - d[3]) / 2e-4;
G[1] = (d[2] - d[4]) / 2e-4;
G[3] = (d[1] + d[3] - 2*currentW) / 1e-8;
G[5] = (d[0] - d[1] - d[2] + 2*currentW - d[3] - d[4] + d[5]) / 2e-8;
G[4] = (d[1] + d[3] - 2*currentW) / 1e-8;
G[2] = (d[1] + d[3] - 2*currentW) / 1e-8;
G[3] = (d[0] - d[1] - d[2] + 2*currentW - d[3] - d[4] + d[5]) / 2e-8;
G[4] = (d[2] + d[4] - 2*currentW) / 1e-8;
// Gradient already zero
if (G[0]*G[0] + G[1]*G[1] > gradTol())
{
//needs to optimise
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
Array<OneD, NekDouble> sk(2), dk(2);
bool DNC = false;
......@@ -177,6 +192,7 @@ void NodeOpti2D2D::Optimise()
// Update node
node->m_x = xc + alpha * sk[0];
node->m_y = yc + alpha * sk[1];
ReCalcMaps();
newVal = GetFunctional<2>(false,false);
//dont need the hessian again this function updates G to be the new
......@@ -207,6 +223,7 @@ void NodeOpti2D2D::Optimise()
//choose whether to do forward or reverse line search
node->m_x = xc + dk[0];
node->m_y = yc + dk[1];
ReCalcMaps();
newVal = GetFunctional<2>(false,false);
if(newVal <= currentW + c1() * (
......@@ -218,11 +235,13 @@ void NodeOpti2D2D::Optimise()
// Update node
node->m_x = xc + alpha * dk[0];
node->m_y = yc + alpha * dk[1];
ReCalcMaps();
newVal = GetFunctional<2>(false,false);
node->m_x = xc + alpha/beta * dk[0];
node->m_y = yc + alpha/beta * dk[1];
ReCalcMaps();
NekDouble dbVal = GetFunctional<2>(false,false);
......@@ -247,6 +266,7 @@ void NodeOpti2D2D::Optimise()
// Update node
node->m_x = xc + alpha * dk[0];
node->m_y = yc + alpha * dk[1];
ReCalcMaps();
newVal = GetFunctional<2>(false,false);
......@@ -267,6 +287,7 @@ void NodeOpti2D2D::Optimise()
//reset the node
node->m_x = xc;
node->m_y = yc;
ReCalcMaps();
mtx.lock();
res->nReset++;
......
......@@ -65,6 +65,8 @@ public:
data[e[i]->GetEl()->GetShapeType()].push_back(e[i]);
}
vertex = false;
//std::map<LibUtilities::ShapeType,std::vector<ElUtilSharedPtr> >::iterator typeIt;
//for(typeIt = data.begin(); typeIt != data.end(); typeIt++)
//{
......@@ -80,6 +82,16 @@ public:
virtual void Optimise() = 0;
NodeOptiJob *GetJob();
void SetVertex()
{
vertex = true;
}
bool IsVertex()
{
return vertex;
}
protected:
template<int DIM> NekDouble GetFunctional(bool gradient = true,
......@@ -89,8 +101,11 @@ protected:
std::map<LibUtilities::ShapeType,std::vector<ElUtilSharedPtr> > data;
Array<OneD, NekDouble> G;
bool vertex;
void CalcMinJac();
bool Linear();
void ReCalcMaps();
template<int DIM> int IsIndefinite();
template<int DIM> void MinEigen(NekDouble &val, Array<OneD, NekDouble> &vec);
......
......@@ -96,15 +96,17 @@ void ProcessVarOpti::BuildDerivUtil()
LibUtilities::ShapeType st = LibUtilities::eQuadrilateral;
derivUtil[st] = boost::shared_ptr<DerivUtil>(new DerivUtil);
derivUtil[st]->ptsLow = m_mesh->m_nummode*m_mesh->m_nummode;
derivUtil[st]->ptsHigh = (m_mesh->m_nummode+4)*(m_mesh->m_nummode+4);
LibUtilities::PointsKey pkey1(m_mesh->m_nummode,
LibUtilities::eGaussLobattoLegendre);
LibUtilities::PointsKey pkey2(m_mesh->m_nummode + 4,
LibUtilities::eNodalQuadSPI);
Array<OneD, NekDouble> u1(derivUtil[st]->ptsLow), v1(derivUtil[st]->ptsLow), u2, v2, tmp;
LibUtilities::eGaussLobattoLegendre);
Array<OneD, NekDouble> u1(derivUtil[st]->ptsLow), v1(derivUtil[st]->ptsLow),
u2(derivUtil[st]->ptsHigh), v2(derivUtil[st]->ptsHigh), tmp, tmp2;
LibUtilities::PointsManager()[pkey1]->GetPoints(tmp);
LibUtilities::PointsManager()[pkey2]->GetPoints(u2, v2);
LibUtilities::PointsManager()[pkey2]->GetPoints(tmp2);
//Get curved nodes ordering for a Quadrilateral is left to right
//ordering so the we need to use gll and redistrobute it.
......@@ -117,7 +119,25 @@ void ProcessVarOpti::BuildDerivUtil()
}
}
derivUtil[st]->ptsHigh = u2.num_elements();
for(int i = 0, ct = 0; i < m_mesh->m_nummode+4; i++)
{
for(int j = 0; j < m_mesh->m_nummode+4; j++, ct++)
{
u2[ct] = tmp2[j];
v2[ct] = tmp2[i];
}
}
Array<OneD, NekDouble> qtmp = LibUtilities::PointsManager()[pkey2]->GetW();
Array<OneD, NekDouble> qds(derivUtil[st]->ptsHigh);
for(int i = 0, ct = 0; i < m_mesh->m_nummode+4; i++)
{
for(int j = 0; j < m_mesh->m_nummode+4; j++, ct++)
{
qds[ct] = qtmp[i]*qtmp[j];
cout << qds[ct] << endl;
}
}
LibUtilities::NodalUtilQuad nodalQuad(
m_mesh->m_nummode - 1, u1, v1);
......@@ -139,7 +159,7 @@ void ProcessVarOpti::BuildDerivUtil()
derivUtil[st]->VdmD[0] = interp * derivUtil[st]->VdmDL[0];
derivUtil[st]->VdmD[1] = interp * derivUtil[st]->VdmDL[1];
//derivUtil->quadW = LibUtilities::MakeQuadratureWeights(U2,V1);
Array<OneD, NekDouble> qds = LibUtilities::PointsManager()[pkey2]->GetW();
//Array<OneD, NekDouble> qds = LibUtilities::PointsManager()[pkey2]->GetW();
NekVector<NekDouble> quadWi(qds);
derivUtil[st]->quadW = quadWi;
......
......@@ -208,6 +208,11 @@ void ProcessVarOpti::Process()
ns.push_back(
GetNodeOptiFactory().CreateInstance(
optiType, freenodes[i][j], it->second, res, derivUtil, opti));
if(freenodes[i][j]->m_id < m_mesh->m_vertexSet.size())
{
ns.back()->SetVertex();
}
}
optiNodes.push_back(ns);
}
......
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