49 namespace NekMeshUtils
54 map<pair<int, int>,
int> nodeorder;
55 map<int, pair<int, int> > nodeorderRev;
73 for (
int i = 0; i < nq - 2; i++)
77 nodeorder[id] = i + 3;
78 nodeorderRev[i + 3] = id;
80 for (
int i = 0; i < nq - 2; i++)
82 id.first = nq - 2 - i;
84 nodeorder[id] = nq + 1 + i;
85 nodeorderRev[nq + 1 + i] = id;
87 for (
int i = 0; i < nq - 2; i++)
90 id.second = nq - 2 - i;
91 nodeorder[id] = nq + nq - 1 + i;
92 nodeorderRev[nq + nq - 1 + i] = id;
98 for (
int k = 0; k < (nq - 3) * (nq - 2) / 2; k++)
102 nodeorder[id] = 3 * (nq - 1) + k;
103 nodeorderRev[3 * (nq - 1) + k] =
id;
113 map<int, vector<int> > nodetosix;
115 for (
int i = (nq + 1) * nq / 2 - (nq - 3) * (nq - 2) / 2;
116 i < (nq + 1) * nq / 2;
122 pair<int, int> p = nodeorderRev[i];
163 set<pair<int, int> > ret;
165 for (it = nodetosix.begin(); it != nodetosix.end(); it++)
167 vector<int> ns = it->second;
168 for (
int i = 0; i < ns.size(); i++)
170 pair<int, int> sp(min(it->first, ns[i]), max(it->first, ns[i]));
192 A(1, 2) = -1.0 + sqrt(3.0);
206 DNekMat C(3, u.num_elements(), 1.0);
207 for (
int i = 0; i < u.num_elements(); i++)
222 for (it = springs.begin(); it != springs.end(); it++)
224 ret[(*it)] = sqrt((pts(0, (*it).first) - pts(0, (*it).second)) *
225 (pts(0, (*it).first) - pts(0, (*it).second)) +
226 (pts(1, (*it).first) - pts(1, (*it).second)) *
227 (pts(1, (*it).first) - pts(1, (*it).second)));
229 if ((*it).first == 12 && (*it).second == 13)
232 if ((*it).first == 12 && (*it).second == 14)
235 if ((*it).first == 13 && (*it).second == 14)
243 if (m_mesh->m_verbose)
244 cout << endl <<
"\tHigh-Order Surface meshing" << endl;
259 int nq = m_mesh->m_nummode;
272 for (
int i = 0; i < m_mesh->m_element[2].size(); i++)
274 if (m_mesh->m_verbose)
277 i, m_mesh->m_element[2].size(),
"\t\tSurface elements");
280 if (m_mesh->m_element[2][i]->GetConf().m_e ==
287 int surf = m_mesh->m_element[2][i]->CADSurfId;
291 vector<EdgeSharedPtr> surfedges =
292 m_mesh->m_element[2][i]->GetEdgeList();
294 vector<EdgeSharedPtr> edges = f->m_edgeList;
295 for (
int j = 0; j < edges.size(); j++)
303 if (test != completedEdges.end())
315 bool foundsurfaceedge =
false;
316 for (
int k = 0; k < surfedges.size(); k++)
318 if (surfedges[k] == e)
320 e->onCurve = surfedges[k]->onCurve;
321 e->CADCurveId = surfedges[k]->CADCurveId;
322 e->CADCurve = surfedges[k]->CADCurve;
323 foundsurfaceedge =
true;
327 "cannot find corresponding surface edge");
331 int cid = e->CADCurveId;
333 NekDouble tb = e->m_n1->GetCADCurveInfo(cid);
334 NekDouble te = e->m_n2->GetCADCurveInfo(cid);
338 for (
int k = 0; k < m_mesh->m_nummode; k++)
341 tb * (1.0 - gll[k]) / 2.0 + te * (1.0 + gll[k]) / 2.0;
345 for (
int k = 1; k < nq - 1; k++)
354 nq - 2, nq - 2, 0.0);
355 for (
int k = 0; k < nq - 2; k++)
362 for (
int k = 0; k < nq - 2; k++)
376 for (
int k = 0; k < nq - 2; k++)
378 Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
390 cout <<
"failed to optimise on curve" << endl;
391 for (
int k = 0; k < nq; k++)
393 ti[k] = tb * (1.0 - gll[k]) / 2.0 +
394 te * (1.0 + gll[k]) / 2.0;
403 cout <<
"BFGS reported no update, curve on "
404 << c->GetId() << endl;
409 ti = opti->GetSolution();
411 vector<CADSurfSharedPtr> s = c->GetAdjSurf();
413 ASSERTL0(s.size() == 2,
"Number of common surfs should be 2");
415 vector<NodeSharedPtr> honodes(m_mesh->m_nummode - 2);
416 for (
int k = 1; k < m_mesh->m_nummode - 1; k++)
420 new Node(0, loc[0], loc[1], loc[2]));
422 nn->SetCADCurve(cid, c, ti[k]);
424 nn->SetCADSurf(s[0]->GetId(), s[0], uv);
425 uv = s[1]->locuv(loc);
426 nn->SetCADSurf(s[1]->GetId(), s[1], uv);
430 e->m_edgeNodes = honodes;
432 completedEdges.insert(e);
439 uvb = e->m_n1->GetCADSurfInfo(surf);
440 uve = e->m_n2->GetCADSurfInfo(surf);
442 for (
int k = 0; k < nq; k++)
445 uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
446 uve[0] * (1.0 + gll[k]) / 2.0;
447 uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
448 uve[1] * (1.0 + gll[k]) / 2.0;
454 for (
int k = 0; k < nq; k++)
456 all[k * 2 + 0] = uvi[k][0];
457 all[k * 2 + 1] = uvi[k][1];
461 for (
int k = 1; k < nq - 1; k++)
463 xi[(k - 1) * 2 + 0] = all[k * 2 + 0];
464 xi[(k - 1) * 2 + 1] = all[k * 2 + 1];
473 for (
int k = 0; k < 2 * (nq - 2); k++)
480 for (
int k = 0; k < 2 * (nq - 2); k++)
492 for (
int k = 0; k < 2 * (nq - 2); k++)
496 Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
501 Norm += J(k, 0) * J(k, 0) / (bnds[3] - bnds[2]) /
515 cout <<
"failed to optimise on edge" << endl;
516 for (
int k = 0; k < nq; k++)
519 uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
520 uve[0] * (1.0 + gll[k]) / 2.0;
521 uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
522 uve[1] * (1.0 + gll[k]) / 2.0;
532 cout <<
"BFGS reported no update, edge on " << surf
539 all = opti->GetSolution();
542 for (
int k = 0; k < nq; k++)
544 uvi[k][0] = all[k * 2 + 0];
545 uvi[k][1] = all[k * 2 + 1];
548 vector<NodeSharedPtr> honodes(nq - 2);
549 for (
int k = 1; k < nq - 1; k++)
554 new Node(0, loc[0], loc[1], loc[2]));
555 nn->SetCADSurf(s->GetId(), s, uvi[k]);
559 e->m_edgeNodes = honodes;
561 completedEdges.insert(e);
565 ASSERTL0(nq <= 5,
"not setup for high-orders yet");
688 if (m_mesh->m_verbose)
#define ASSERTL0(condition, msg)
boost::shared_ptr< OptiEdge > OptiEdgeSharedPtr
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Represents a point in the domain.
map< pair< int, int >, NekDouble > weights(set< pair< int, int > > springs, Array< OneD, NekDouble > u, Array< OneD, NekDouble > v)
set< pair< int, int > > ListOfFaceSpings(int nq)
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
2D Nodal Fekete Points on a Triangle
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
boost::shared_ptr< CADSurf > CADSurfSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
bool BGFSUpdate(OptiObjSharedPtr opti, DNekMat &J, DNekMat &B, DNekMat &H)
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
HOTriangle< NodeSharedPtr > HOSurf
void PrintProgressbar(const int position, const int goal, const string message)
Prints a progressbar.
1D Gauss-Lobatto-Legendre quadrature points
boost::shared_ptr< CADCurve > CADCurveSharedPtr