Commit 8f29693f authored by Michael Turner's avatar Michael Turner
Browse files

bug in curve optimiser, could be due to DNC applicability

parent 03563d1c
......@@ -214,9 +214,9 @@ void Mesh::MakeOrder(int order,
}
// Copy face curvature
el->SetVolumeNodes(edge->m_edgeNodes);
el->MakeOrder(order, SpatialDomains::GeometrySharedPtr(),
pTypes[el->GetConf().m_e], m_spaceDim, id, true);
el->SetVolumeNodes(edge->m_edgeNodes);
}
for (int i = 0; i < m_element[2].size(); ++i)
......
......@@ -86,7 +86,7 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
if(m_el->GetConf().m_e == LibUtilities::eQuadrilateral)
{
vector<Array<OneD, NekDouble> > xyz(6);
vector<Array<OneD, NekDouble> > xyz(4);
vector<NodeSharedPtr> ns = m_el->GetVertexList();
for(int i = 0; i < 4; i++)
{
......@@ -97,8 +97,6 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
xyz[i] = x;
}
vector<DNekMat> tmp;
for (int i = 0; i < derivUtil->ptsHigh; ++i)
{
NekDouble a1 = 0.5 * (1 - derivUtil->ptx[i]);
......@@ -120,8 +118,6 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
J.Invert();
tmp.push_back(J);
Array<OneD, NekDouble> r(10,0.0); //store det in 10th entry
r[9] = 1.0 / (J(0,0) * J(1,1) - J(0,1) * J(1,0));
......@@ -226,26 +222,6 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
xyz[i] = x;
}
ElmtConfig c = m_el->GetConf();
c.m_order = 1;
c.m_reorient = false;
c.m_volumeNodes = false;
c.m_faceNodes = false;
ElementSharedPtr E = GetElementFactory().CreateInstance(
c.m_e, c, ns, m_el->GetTagList());
SpatialDomains::GeometrySharedPtr geom = m_el->GetGeom(3);
LibUtilities::PointsKeyVector p = geom->GetPointsKeys();
SpatialDomains::GeomFactorsSharedPtr gfac = geom->GetGeomFactors();
Array<OneD, NekDouble> jc = gfac->GetJac(p);
StdRegions::StdExpansionSharedPtr xmap = geom->GetXmap();
Array<OneD, NekDouble> coeff(xmap->GetNcoeffs());
xmap->FwdTrans(jc,coeff);
Array<OneD, NekDouble> phys(xmap->GetTotPoints());
xmap->BwdTrans(coeff, phys);
for (int i = 0; i < derivUtil->ptsHigh; ++i)
{
NekDouble a2 = 0.5 * (1 + derivUtil->ptx[i]);
......@@ -288,13 +264,6 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
-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)));
Array<OneD, NekDouble> xp(3);
xp[0] = derivUtil->ptx[i];
xp[1] = derivUtil->pty[i];
xp[2] = derivUtil->ptz[i];
cout << r[9] << " " << xmap->PhysEvaluate(xp, phys) << endl;
r[0] = J(0,0);
r[1] = J(1,0);
r[2] = J(2,0);
......@@ -306,7 +275,6 @@ vector<Array<OneD, NekDouble> > ElUtil::MappingIdealToRef()
r[8] = J(2,2);
ret.push_back(r);
}
exit(-1);
}
else if (m_el->GetConf().m_e == LibUtilities::eHexahedron)
......@@ -446,21 +414,19 @@ void ElUtil::Evaluate()
Z(j) = *nodes[j][2];
}
NekVector<NekDouble> x1i(nodes.size()),y1i(nodes.size()),z1i(nodes.size()),
x2i(nodes.size()),y2i(nodes.size()),z2i(nodes.size()),
x3i(nodes.size()),y3i(nodes.size()),z3i(nodes.size());
x1i = derivUtil->VdmDL[0]*X;
y1i = derivUtil->VdmDL[0]*Y;
z1i = derivUtil->VdmDL[0]*Z;
x2i = derivUtil->VdmDL[1]*X;
y2i = derivUtil->VdmDL[1]*Y;
z2i = derivUtil->VdmDL[1]*Z;
x3i = derivUtil->VdmDL[2]*X;
y3i = derivUtil->VdmDL[2]*Y;
z3i = derivUtil->VdmDL[2]*Z;
NekVector<NekDouble> x1i(derivUtil->ptsHigh),y1i(derivUtil->ptsHigh),z1i(derivUtil->ptsHigh),
x2i(derivUtil->ptsHigh),y2i(derivUtil->ptsHigh),z2i(derivUtil->ptsHigh),
x3i(derivUtil->ptsHigh),y3i(derivUtil->ptsHigh),z3i(derivUtil->ptsHigh);
Array<OneD, NekDouble> jacs(nodes.size());
x1i = derivUtil->VdmD[0]*X;
y1i = derivUtil->VdmD[0]*Y;
z1i = derivUtil->VdmD[0]*Z;
x2i = derivUtil->VdmD[1]*X;
y2i = derivUtil->VdmD[1]*Y;
z2i = derivUtil->VdmD[1]*Z;
x3i = derivUtil->VdmD[2]*X;
y3i = derivUtil->VdmD[2]*Y;
z3i = derivUtil->VdmD[2]*Z;
/*if(m_el->GetShapeType() == LibUtilities::ePrism)
{
......@@ -497,7 +463,7 @@ void ElUtil::Evaluate()
}
}*/
for(int j = 0; j < nodes.size(); j++)
for(int j = 0; j < derivUtil->ptsHigh; j++)
{
DNekMat dxdz(3,3,1.0,eFULL);
dxdz(0,0) = x1i(j);
......
......@@ -191,6 +191,8 @@ inline NekDouble CalcIdealJac(int elmt,
}
}
//this matrix multiplier has been checked MT 14/9
return Determinant(jacIdeal);
}
......@@ -258,7 +260,7 @@ NekDouble NodeOpti::GetFunctional(bool gradient, bool hessian)
NekDouble integral = 0.0;
NekDouble ep = minJac < 1e-7 ? sqrt(1e-14 + 0.004*minJac*minJac) : 1e-7;
NekDouble ep = minJac < 1e-4 ? sqrt(1e-8 + 0.004*minJac*minJac) : 1e-4;
//NekDouble ep = minJac < 0.0 ? sqrt(gam*(gam-minJac)) : gam;
//NekDouble ep = minJac < gam ? sqrt(gam*gam + minJac*minJac) : gam;
NekDouble jacIdeal[DIM][DIM], jacDet;
......@@ -430,15 +432,25 @@ NekDouble NodeOpti::GetFunctional(bool gradient, bool hessian)
jacDet = CalcIdealJac(i, k, derivs[typeIt->first], typeIt->second, jacIdeal);
NekDouble I1 = FrobeniusNorm(jacIdeal);
//if(gradient && minJac > 0.0 && jacDet < 0.0)
//{
// ASSERTL0(false,"minjac is positive but jacdet is negative");
// return numeric_limits<double>::max();
//}
//ep = minJac < 0.0 ? sqrt(gam*jacDet*(gam*jacDet-minJac)) : gam*jacDet;
NekDouble sigma =
0.5*(jacDet + sqrt(jacDet*jacDet + 4.0*ep*ep));
if(sigma < 1e-20)
if(sigma < numeric_limits<float>::min())
{
sigma = 1e-20;
sigma = numeric_limits<float>::min();
}
NekDouble lsigma = log(sigma);
integral += derivUtil[typeIt->first]->quadW[k]*
fabs(typeIt->second[i]->maps[k][9]) *
fabs(typeIt->second[i]->maps[k][9]) *
(0.5 * mu * (I1 - 3.0 - 2.0*lsigma) +
0.5 * K * lsigma * lsigma);
......
......@@ -95,15 +95,36 @@ void NodeOpti2D2D::Optimise()
{
CalcMinJac();
NekDouble currentW = GetFunctional<2>();
NekDouble currentW = GetFunctional<2>(false,false);
NekDouble newVal = currentW;
return;
NekDouble xc = node->m_x;
NekDouble yc = node->m_y;
vector<NekDouble> d;
for(int i = 0; i < 6; i++)
{
node->m_x = xc + 1e-6*dir[i][0];
node->m_y = yc + 1e-6*dir[i][1];
d.push_back(GetFunctional<2>(false,false));
}
G = Array<OneD, NekDouble>(5);
G[0] = (d[1] - d[3]) / 2e-6;
G[1] = (d[2] - d[4]) / 2e-6;
G[2] = (d[1] + d[3] - 2*currentW) / 1e-12;
G[3] = (d[0] - d[1] - d[2] + 2*currentW - d[3] - d[4] + d[5]) / 2e-12;
G[4] = (d[2] + d[4] - 2*currentW) / 1e-12;
// 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;
//NekDouble xc = node->m_x;
//NekDouble yc = node->m_y;
Array<OneD, NekDouble> sk(2), dk(2);
bool DNC = false;
......@@ -462,9 +483,11 @@ void NodeOpti3D3D::Optimise()
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);
res->func +=newVal;
mtx.unlock();
}
mtx.lock();
res->func += newVal;
mtx.unlock();
}
NodeOptiJob* NodeOpti::GetJob()
......
......@@ -56,11 +56,15 @@ void NodeOpti1D3D::Optimise()
NekDouble currentW = GetFunctional<3>();
NekDouble newVal = currentW;
cout << "curve " << sqrt(G[0]*G[0] + G[1]*G[1] + G[2]*G[2]) << endl;
if (G[0]*G[0] + G[1]*G[1] + G[2]*G[2] > gradTol())
{
//modify the gradient to be on the cad system
ProcessGradient();
cout << G[0] << " " << G[1] << endl;
//needs to optimise
NekDouble tc = node->GetCADCurveInfo(curve->GetId());
NekDouble xc = node->m_x;
......@@ -277,13 +281,16 @@ void NodeOpti1D3D::Optimise()
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);
res->func += newVal;
mtx.unlock();
}
mtx.lock();
res->func += newVal;
mtx.unlock();
}
int NodeOpti2D3D::m_type = GetNodeOptiFactory().RegisterCreatorFunction(
23, NodeOpti2D3D::create, "1D3D");
23, NodeOpti2D3D::create, "2D3D");
void NodeOpti2D3D::Optimise()
{
......@@ -292,6 +299,8 @@ void NodeOpti2D3D::Optimise()
NekDouble currentW = GetFunctional<3>();
NekDouble newVal = currentW;
//cout << "curve " << sqrt(G[0]*G[0] + G[1]*G[1] + G[2]*G[2]) << endl;
if (G[0]*G[0] + G[1]*G[1] + G[2]*G[2] > gradTol())
{
//modify the gradient to be on the cad system
......@@ -529,7 +538,6 @@ void NodeOpti2D3D::Optimise()
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);
res->func += newVal;
mtx.unlock();
Array<OneD, NekDouble> uva = node->GetCADSurfInfo(surf->GetId());
......@@ -539,6 +547,10 @@ void NodeOpti2D3D::Optimise()
ASSERTL1(false,"something when very wrong and the node finished out of its parameter plane");
}
}
mtx.lock();
res->func += newVal;
mtx.unlock();
}
void NodeOpti1D3D::ProcessGradient()
......
......@@ -694,7 +694,7 @@ vector<ElementSharedPtr> ProcessVarOpti::GetLockedElements(NekDouble thres)
}
}
for (int i = 0; i < 12; i++)
for (int i = 0; i < 8; i++)
{
vector<ElementSharedPtr> tmp = totest;
totest.clear();
......
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