Commit 0d4ed54b authored by Michael Turner's avatar Michael Turner
Browse files

fix surface curve bug

parent 8c18fdf9
......@@ -126,7 +126,7 @@ public:
/// apply spherigon surface smoothing.
std::set<std::pair<int, int> > m_spherigonSurfs;
/// List of face labels for composite annotation
std::map<int, std::string> m_faceLabels;
std::map<int, std::string> m_faceLabels;
/// Whether the mesh has CAD
bool m_hasCAD;
/// CAD file ID
......
......@@ -327,7 +327,7 @@ void Module::ProcessFaces(bool ReprocessFaces)
for (int j = 0; j < elmt->GetVertexCount(); ++j)
{
elmt->SetVertex(j, (*it)->m_vertexList[j], false);
elmt->SetEdge(j, (*it)->m_edgeList[j], false);
//elmt->SetEdge(j, (*it)->m_edgeList[j], false);
}
// Update 3D element boundary map.
......
......@@ -213,7 +213,7 @@ void CurveMesh::Mesh()
new Edge(m_meshpoints[i], m_meshpoints[i + 1]));
e->CADCurveId = m_id;
e->CADCurve = m_cadcurve;
e->onCurve = true;
e->onCurve = true;
m_mesh->m_edgeSet.insert(e);
}
......
......@@ -828,7 +828,7 @@ void FaceMesh::BuildLocalMesh()
if (!(s == m_mesh->m_edgeSet.end()))
{
edgs[j] = *s;
E->SetEdge(j, edgs[j]);
E->SetEdge(j, *s);
}
pair<EdgeSet::iterator, bool> test = m_localEdges.insert(edgs[j]);
......
......@@ -215,12 +215,6 @@ inline map<pair<int, int>, NekDouble> weights(set<pair<int, int> > springs,
DNekMat pts = M * C;
/*for(int i = 0; i < u.num_elements(); i++)
{
cout << pts(0,i) << " " << pts(1,i) << endl;
}
exit(-1);*/
set<pair<int, int> >::iterator it;
for (it = springs.begin(); it != springs.end(); it++)
{
......@@ -228,15 +222,6 @@ inline map<pair<int, int>, NekDouble> weights(set<pair<int, int> > springs,
(pts(0, (*it).first) - pts(0, (*it).second)) +
(pts(1, (*it).first) - pts(1, (*it).second)) *
(pts(1, (*it).first) - pts(1, (*it).second)));
if ((*it).first == 12 && (*it).second == 13)
ret[(*it)] *= 1.2;
if ((*it).first == 12 && (*it).second == 14)
ret[(*it)] *= 1.2;
if ((*it).first == 13 && (*it).second == 14)
ret[(*it)] *= 1.2;
}
return ret;
}
......@@ -248,8 +233,8 @@ ModuleKey HOSurfaceMesh::className = GetModuleFactory().RegisterCreatorFunction(
HOSurfaceMesh::HOSurfaceMesh(MeshSharedPtr m) : ProcessModule(m)
{
m_config["extract"] =
ConfigOption(true, "0", "Extract non-valid elements from mesh.");
m_config["opti"] =
ConfigOption(true, "0", "Perform edge node optimisation.");
}
HOSurfaceMesh::~HOSurfaceMesh()
......@@ -290,12 +275,21 @@ void HOSurfaceMesh::Process()
set<pair<int, int> > springs = ListOfFaceSpings(nq);
map<pair<int, int>, NekDouble> z = weights(springs, u, v);
// because edges are listed twice need a method to not repeat over them
EdgeSet completedEdges;
bool qOpti = m_config["opti"].beenSet;
// loop over all the faces in the surface mesh, check all three edges for
// high order info, if nothing high-order the edge.
EdgeSet surfaceEdges;
EdgeSet completedEdges;
for(int i = 0; i < m_mesh->m_element[2].size(); i++)
{
vector<EdgeSharedPtr> es = m_mesh->m_element[2][i]->GetEdgeList();
for(int j = 0; j < es.size(); j++)
surfaceEdges.insert(es[j]);
}
for (int i = 0; i < m_mesh->m_element[2].size(); i++)
{
if (m_mesh->m_verbose)
......@@ -313,43 +307,33 @@ void HOSurfaceMesh::Process()
f->CADSurfId = surf;
f->CADSurf = s;
vector<EdgeSharedPtr> surfedges =
m_mesh->m_element[2][i]->GetEdgeList();
vector<EdgeSharedPtr> edges = f->m_edgeList;
for (int j = 0; j < edges.size(); j++)
{
EdgeSharedPtr e = edges[j];
// test insert the edge into completedEdges
// if the edge already exists move on
// if not figure out its high-order information
EdgeSet::iterator test = completedEdges.find(edges[j]);
EdgeSet::iterator test = completedEdges.find(e);
if (test != completedEdges.end())
{
continue;
}
EdgeSharedPtr e = edges[j];
// the edges in the element are different to those in the face
// the cad information is stored in the element edges which are not
// in the m_mesh->m_edgeSet group.
// in the m_mesh->m_edgeSet groups
// need to link them together and copy the cad information to be
// able to identify how to make it high-order
bool foundsurfaceedge = false;
for (int k = 0; k < surfedges.size(); k++)
{
if (surfedges[k] == e)
{
e->onCurve = surfedges[k]->onCurve;
e->CADCurveId = surfedges[k]->CADCurveId;
e->CADCurve = surfedges[k]->CADCurve;
foundsurfaceedge = true;
}
}
ASSERTL0(foundsurfaceedge,
"cannot find corresponding surface edge");
EdgeSet::iterator it = surfaceEdges.find(e);
ASSERTL0(it != surfaceEdges.end(),"could not find edge in surface");
e->onCurve = (*it)->onCurve;
e->CADCurveId = (*it)->CADCurveId;
e->CADCurve = (*it)->CADCurve;
vector<NodeSharedPtr> honodes(m_mesh->m_nummode - 2);
if (e->onCurve)
{
......@@ -366,80 +350,82 @@ void HOSurfaceMesh::Process()
tb * (1.0 - gll[k]) / 2.0 + te * (1.0 + gll[k]) / 2.0;
}
/*Array<OneD, NekDouble> xi(nq - 2);
for (int k = 1; k < nq - 1; k++)
{
xi[k - 1] = ti[k];
}
OptiEdgeSharedPtr opti =
MemoryManager<OptiEdge>::AllocateSharedPtr(ti, gll, c);
DNekMat B(
nq - 2, nq - 2, 0.0); // approximate hessian (I to start)
for (int k = 0; k < nq - 2; k++)
{
B(k, k) = 1.0;
}
DNekMat H(nq - 2,
nq - 2,
0.0); // approximate inverse hessian (I to start)
for (int k = 0; k < nq - 2; k++)
if(qOpti)
{
H(k, k) = 1.0;
}
DNekMat J = opti->dF(xi);
Array<OneD, NekDouble> xi(nq - 2);
for (int k = 1; k < nq - 1; k++)
{
xi[k - 1] = ti[k];
}
Array<OneD, NekDouble> bnds = c->Bounds();
OptiEdgeSharedPtr opti =
MemoryManager<OptiEdge>::AllocateSharedPtr(ti, gll, c);
bool repeat = true;
int itct = 0;
while (repeat)
{
NekDouble Norm = 0;
DNekMat B(
nq - 2, nq - 2, 0.0); // approximate hessian (I to start)
for (int k = 0; k < nq - 2; k++)
{
Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
(bnds[1] - bnds[0]);
B(k, k) = 1.0;
}
Norm = sqrt(Norm);
if (Norm < 1E-8)
DNekMat H(nq - 2,
nq - 2,
0.0); // approximate inverse hessian (I to start)
for (int k = 0; k < nq - 2; k++)
{
repeat = false;
break;
H(k, k) = 1.0;
}
if (itct > 1000)
DNekMat J = opti->dF(xi);
Array<OneD, NekDouble> bnds = c->Bounds();
bool repeat = true;
int itct = 0;
while (repeat)
{
cout << "failed to optimise on curve" << endl;
for (int k = 0; k < nq; k++)
NekDouble Norm = 0;
for (int k = 0; k < nq - 2; k++)
{
ti[k] = tb * (1.0 - gll[k]) / 2.0 +
te * (1.0 + gll[k]) / 2.0;
Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
(bnds[1] - bnds[0]);
}
break;
}
itct++;
Norm = sqrt(Norm);
if (!BGFSUpdate(opti, J, B, H))
{
if(m_mesh->m_verbose)
if (Norm < 1E-8)
{
repeat = false;
break;
}
if (itct > 1000)
{
cout << "BFGS reported no update, curve on "
<< c->GetId() << endl;
cout << "failed to optimise on curve" << endl;
for (int k = 0; k < nq; k++)
{
ti[k] = tb * (1.0 - gll[k]) / 2.0 +
te * (1.0 + gll[k]) / 2.0;
}
break;
}
itct++;
if (!BGFSUpdate(opti, J, B, H))
{
if(m_mesh->m_verbose)
{
cout << "BFGS reported no update, curve on "
<< c->GetId() << endl;
}
break;
}
break;
}
// need to pull the solution out of opti
ti = opti->GetSolution();
}
// need to pull the solution out of opti
ti = opti->GetSolution();*/
vector<CADSurfSharedPtr> s = c->GetAdjSurf();
ASSERTL0(s.size() == 2, "Number of common surfs should be 2");
vector<NodeSharedPtr> honodes(m_mesh->m_nummode - 2);
for (int k = 1; k < m_mesh->m_nummode - 1; k++)
{
Array<OneD, NekDouble> loc = c->P(ti[k]);
......@@ -453,10 +439,6 @@ void HOSurfaceMesh::Process()
nn->SetCADSurf(s[1]->GetId(), s[1], uv);
honodes[k - 1] = nn;
}
e->m_edgeNodes = honodes;
e->m_curveType = LibUtilities::eGaussLobattoLegendre;
completedEdges.insert(e);
}
else
{
......@@ -479,105 +461,107 @@ void HOSurfaceMesh::Process()
uvi[k] = uv;
}
/*Array<OneD, NekDouble> bnds = s->GetBounds();
Array<OneD, NekDouble> all(2 * nq);
for (int k = 0; k < nq; k++)
if(qOpti)
{
all[k * 2 + 0] = uvi[k][0];
all[k * 2 + 1] = uvi[k][1];
}
Array<OneD, NekDouble> bnds = s->GetBounds();
Array<OneD, NekDouble> all(2 * nq);
for (int k = 0; k < nq; k++)
{
all[k * 2 + 0] = uvi[k][0];
all[k * 2 + 1] = uvi[k][1];
}
Array<OneD, NekDouble> xi(2 * (nq - 2));
for (int k = 1; k < nq - 1; k++)
{
xi[(k - 1) * 2 + 0] = all[k * 2 + 0];
xi[(k - 1) * 2 + 1] = all[k * 2 + 1];
}
Array<OneD, NekDouble> xi(2 * (nq - 2));
for (int k = 1; k < nq - 1; k++)
{
xi[(k - 1) * 2 + 0] = all[k * 2 + 0];
xi[(k - 1) * 2 + 1] = all[k * 2 + 1];
}
OptiEdgeSharedPtr opti =
MemoryManager<OptiEdge>::AllocateSharedPtr(all, gll, s);
OptiEdgeSharedPtr opti =
MemoryManager<OptiEdge>::AllocateSharedPtr(all, gll, s);
DNekMat B(2 * (nq - 2),
2 * (nq - 2),
0.0); // approximate hessian (I to start)
for (int k = 0; k < 2 * (nq - 2); k++)
{
B(k, k) = 1.0;
}
DNekMat H(2 * (nq - 2),
2 * (nq - 2),
0.0); // approximate inverse hessian (I to start)
for (int k = 0; k < 2 * (nq - 2); k++)
{
H(k, k) = 1.0;
}
DNekMat B(2 * (nq - 2),
2 * (nq - 2),
0.0); // approximate hessian (I to start)
for (int k = 0; k < 2 * (nq - 2); k++)
{
B(k, k) = 1.0;
}
DNekMat H(2 * (nq - 2),
2 * (nq - 2),
0.0); // approximate inverse hessian (I to start)
for (int k = 0; k < 2 * (nq - 2); k++)
{
H(k, k) = 1.0;
}
DNekMat J = opti->dF(xi);
DNekMat J = opti->dF(xi);
bool repeat = true;
int itct = 0;
while (repeat)
{
NekDouble Norm = 0;
for (int k = 0; k < 2 * (nq - 2); k++)
bool repeat = true;
int itct = 0;
while (repeat)
{
if (k % 2 == 0)
NekDouble Norm = 0;
for (int k = 0; k < 2 * (nq - 2); k++)
{
Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
(bnds[1] - bnds[0]);
if (k % 2 == 0)
{
Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
(bnds[1] - bnds[0]);
}
else
{
Norm += J(k, 0) * J(k, 0) / (bnds[3] - bnds[2]) /
(bnds[3] - bnds[2]);
}
}
else
Norm = sqrt(Norm);
if (Norm < 1E-8)
{
Norm += J(k, 0) * J(k, 0) / (bnds[3] - bnds[2]) /
(bnds[3] - bnds[2]);
repeat = false;
break;
}
}
Norm = sqrt(Norm);
if (Norm < 1E-8)
{
repeat = false;
break;
}
if (itct > 1000)
{
cout << "failed to optimise on edge" << endl;
for (int k = 0; k < nq; k++)
if (itct > 1000)
{
Array<OneD, NekDouble> uv(2);
uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
uve[0] * (1.0 + gll[k]) / 2.0;
uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
uve[1] * (1.0 + gll[k]) / 2.0;
uvi[k] = uv;
cout << "failed to optimise on edge" << endl;
for (int k = 0; k < nq; k++)
{
Array<OneD, NekDouble> uv(2);
uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
uve[0] * (1.0 + gll[k]) / 2.0;
uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
uve[1] * (1.0 + gll[k]) / 2.0;
uvi[k] = uv;
}
break;
}
break;
}
itct++;
itct++;
if (!BGFSUpdate(opti, J, B, H))
{
if(m_mesh->m_verbose)
if (!BGFSUpdate(opti, J, B, H))
{
cout << "BFGS reported no update, edge on " << surf
<< endl;
if(m_mesh->m_verbose)
{
cout << "BFGS reported no update, edge on " << surf
<< endl;
}
// exit(-1);
break;
}
// exit(-1);
break;
}
}
all = opti->GetSolution();
all = opti->GetSolution();
// need to put all backinto uv
for (int k = 0; k < nq; k++)
{
uvi[k][0] = all[k * 2 + 0];
uvi[k][1] = all[k * 2 + 1];
}*/
// need to put all backinto uv
for (int k = 0; k < nq; k++)
{
uvi[k][0] = all[k * 2 + 0];
uvi[k][1] = all[k * 2 + 1];
}
}
vector<NodeSharedPtr> honodes(nq - 2);
for (int k = 1; k < nq - 1; k++)
{
Array<OneD, NekDouble> loc;
......@@ -587,14 +571,14 @@ void HOSurfaceMesh::Process()
nn->SetCADSurf(s->GetId(), s, uvi[k]);
honodes[k - 1] = nn;
}
e->m_edgeNodes = honodes;
e->m_curveType = LibUtilities::eGaussLobattoLegendre;
completedEdges.insert(e);
}
e->m_edgeNodes = honodes;
e->m_curveType = LibUtilities::eGaussLobattoLegendre;
completedEdges.insert(e);
}
ASSERTL0(nq <= 5, "not setup for high-orders yet");
//just add the face interior nodes through interp and project
vector<NodeSharedPtr> vertices = f->m_vertexList;
......@@ -686,6 +670,7 @@ void HOSurfaceMesh::Process()
f->m_faceNodes = honodes;
f->m_curveType = LibUtilities::eGaussLobattoLegendre;
}
}
if (m_mesh->m_verbose)
......
////////////////////////////////////////////////////////////////////////////////
//
// File: ProcessJac.h
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Calculate jacobians of elements.
//
////////////////////////////////////////////////////////////////////////////////
#ifndef NEKMESHUTILS_SURFACE_HO
#define NEKMESHUTILS_SURFACE_HO
#include <NekMeshUtils/Module/Module.h>
namespace Nektar
{
namespace NekMeshUtils
{
class HOSurfaceMesh : public ProcessModule
{
public:
/// Creates an instance of this class
static boost::shared_ptr<Module> create(MeshSharedPtr m)
{
return MemoryManager<HOSurfaceMesh>::AllocateSharedPtr(m);
}
static ModuleKey className;
HOSurfaceMesh(MeshSharedPtr m);
virtual ~HOSurfaceMesh();
virtual void Process();
};
}
}
#endif
......@@ -89,7 +89,8 @@ void InputCAD::Process()
pSession->LoadParameter("MaxDelta", m_maxDelta);
pSession->LoadParameter("EPS", m_eps);
pSession->LoadParameter("Order", m_order);
m_CADName = pSession->GetSolverInfo("CADFile");
m_mesh->m_CADId = pSession->GetSolverInfo("CADFile");
m_mesh->m_hasCAD = true;
if (pSession->DefinesSolverInfo("MeshType"))
{
......@@ -113,18 +114,17 @@ void InputCAD::Process()
m_writeoctree = pSession->GetSolverInfo("WriteOctree") == "TRUE";
}