73 cout << endl <<
"\tHigh-Order Surface meshing" << endl;
75 LibUtilities::PointsKey ekey(
m_mesh->m_nummode,
77 Array<OneD, NekDouble> gll;
81 LibUtilities::PointsKey pkey(
m_mesh->m_nummode,
84 Array<OneD, NekDouble> u, v;
86 int nq =
m_mesh->m_nummode;
88 int np = nq * (nq + 1) / 2;
90 int ni = (nq-2)*(nq-3) / 2;
94 int niq = npq - 4 - 4*(nq-2);
98 bool qOpti =
m_config[
"opti"].beenSet;
106 for(
int i = 0; i <
m_mesh->m_element[2].size(); i++)
108 vector<EdgeSharedPtr> es =
m_mesh->m_element[2][i]->GetEdgeList();
109 for(
int j = 0; j < es.size(); j++)
110 surfaceEdges.insert(es[j]);
113 for (
int i = 0; i <
m_mesh->m_element[2].size(); i++)
118 i,
m_mesh->m_element[2].size(),
"\t\tSurface elements");
123 int surf = s->GetId();
129 f = boost::shared_ptr<Face>(
new Face(
m_mesh->m_element[2][i]->GetVertexList(),
130 vector<NodeSharedPtr>(),
131 m_mesh->m_element[2][i]->GetEdgeList(),
137 vector<EdgeSharedPtr> edges = f->m_edgeList;
138 for (
int j = 0; j < edges.size(); j++)
147 if (test != completedEdges.end())
158 ASSERTL0(it != surfaceEdges.end(),
"could not find edge in surface");
160 if((*it)->m_parentCAD)
162 e->m_parentCAD = (*it)->m_parentCAD;
169 vector<NodeSharedPtr> honodes(
m_mesh->m_nummode - 2);
173 int cid = e->m_parentCAD->GetId();
175 NekDouble tb = e->m_n1->GetCADCurveInfo(cid);
176 NekDouble te = e->m_n2->GetCADCurveInfo(cid);
179 Array<OneD, NekDouble> ti(
m_mesh->m_nummode);
180 for (
int k = 0; k <
m_mesh->m_nummode; k++)
183 tb * (1.0 - gll[k]) / 2.0 + te * (1.0 + gll[k]) / 2.0;
188 Array<OneD, NekDouble> xi(nq - 2);
189 for (
int k = 1; k < nq - 1; k++)
198 nq - 2, nq - 2, 0.0);
199 for (
int k = 0; k < nq - 2; k++)
206 for (
int k = 0; k < nq - 2; k++)
213 Array<OneD, NekDouble> bnds = c->GetBounds();
220 for (
int k = 0; k < nq - 2; k++)
222 Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
234 cout <<
"failed to optimise on curve" << endl;
235 for (
int k = 0; k < nq; k++)
237 ti[k] = tb * (1.0 - gll[k]) / 2.0 +
238 te * (1.0 + gll[k]) / 2.0;
248 cout <<
"BFGS reported no update, curve on "
249 << c->GetId() << endl;
255 ti = opti->GetSolution();
257 vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s = c->GetAdjSurf();
259 for (
int k = 1; k <
m_mesh->m_nummode - 1; k++)
261 Array<OneD, NekDouble> loc = c->P(ti[k]);
263 new Node(0, loc[0], loc[1], loc[2]));
265 nn->SetCADCurve(cid, c, ti[k]);
266 for(
int m = 0; m < s.size(); m++)
268 Array<OneD, NekDouble> uv = s[m].first->locuv(loc);
269 nn->SetCADSurf(s[m].first->GetId(), s[m].first, uv);
278 Array<OneD, NekDouble> uvb, uve;
279 uvb = e->m_n1->GetCADSurfInfo(surf);
280 uve = e->m_n2->GetCADSurfInfo(surf);
282 Array<OneD, Array<OneD, NekDouble> > uvi(nq);
283 for (
int k = 0; k < nq; k++)
285 Array<OneD, NekDouble> uv(2);
286 uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
287 uve[0] * (1.0 + gll[k]) / 2.0;
288 uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
289 uve[1] * (1.0 + gll[k]) / 2.0;
295 Array<OneD, NekDouble> bnds = s->GetBounds();
296 Array<OneD, NekDouble> all(2 * nq);
297 for (
int k = 0; k < nq; k++)
299 all[k * 2 + 0] = uvi[k][0];
300 all[k * 2 + 1] = uvi[k][1];
303 Array<OneD, NekDouble> xi(2 * (nq - 2));
304 for (
int k = 1; k < nq - 1; k++)
306 xi[(k - 1) * 2 + 0] = all[k * 2 + 0];
307 xi[(k - 1) * 2 + 1] = all[k * 2 + 1];
316 for (
int k = 0; k < 2 * (nq - 2); k++)
323 for (
int k = 0; k < 2 * (nq - 2); k++)
335 for (
int k = 0; k < 2 * (nq - 2); k++)
339 Norm += J(k, 0) * J(k, 0) / (bnds[1] - bnds[0]) /
344 Norm += J(k, 0) * J(k, 0) / (bnds[3] - bnds[2]) /
358 cout <<
"failed to optimise on edge" << endl;
359 for (
int k = 0; k < nq; k++)
361 Array<OneD, NekDouble> uv(2);
362 uv[0] = uvb[0] * (1.0 - gll[k]) / 2.0 +
363 uve[0] * (1.0 + gll[k]) / 2.0;
364 uv[1] = uvb[1] * (1.0 - gll[k]) / 2.0 +
365 uve[1] * (1.0 + gll[k]) / 2.0;
376 cout <<
"BFGS reported no update, edge on " << surf
384 all = opti->GetSolution();
387 for (
int k = 0; k < nq; k++)
389 uvi[k][0] = all[k * 2 + 0];
390 uvi[k][1] = all[k * 2 + 1];
394 for (
int k = 1; k < nq - 1; k++)
396 Array<OneD, NekDouble> loc;
399 new Node(0, loc[0], loc[1], loc[2]));
400 nn->SetCADSurf(s->GetId(), s, uvi[k]);
405 e->m_edgeNodes = honodes;
407 completedEdges.insert(e);
411 vector<NodeSharedPtr> vertices = f->m_vertexList;
416 Array<OneD, NekDouble> coeffs0 = geom->GetCoeffs(0);
417 Array<OneD, NekDouble> coeffs1 = geom->GetCoeffs(1);
418 Array<OneD, NekDouble> coeffs2 = geom->GetCoeffs(2);
420 Array<OneD, NekDouble> xc(xmap->GetTotPoints());
421 Array<OneD, NekDouble> yc(xmap->GetTotPoints());
422 Array<OneD, NekDouble> zc(xmap->GetTotPoints());
424 xmap->BwdTrans(coeffs0,xc);
425 xmap->BwdTrans(coeffs1,yc);
426 xmap->BwdTrans(coeffs2,zc);
428 if(vertices.size() == 3)
431 vector<Array<OneD, NekDouble> > uvi;
432 for(
int j = np-ni; j < np; j++)
434 Array<OneD, NekDouble> xp(2);
438 Array<OneD, NekDouble> loc(3);
439 loc[0] = xmap->PhysEvaluate(xp, xc);
440 loc[1] = xmap->PhysEvaluate(xp, yc);
441 loc[2] = xmap->PhysEvaluate(xp, zc);
443 Array<OneD, NekDouble> uv(2);
444 s->ProjectTo(loc,uv);
449 vector<NodeSharedPtr> honodes;
450 for(
int j = 0; j < ni; j++)
452 Array<OneD, NekDouble> loc;
455 Node(0,loc[0],loc[1],loc[2]));
456 nn->SetCADSurf(surf, s, uvi[j]);
457 honodes.push_back(nn);
460 f->m_faceNodes = honodes;
463 else if(vertices.size() == 4)
466 vector<Array<OneD, NekDouble> > uvi;
467 for(
int i = 1; i < nq - 1; i++)
469 for(
int j = 1; j < nq - 1; j++)
471 Array<OneD, NekDouble> xp(2);
475 Array<OneD, NekDouble> loc(3);
476 loc[0] = xmap->PhysEvaluate(xp, xc);
477 loc[1] = xmap->PhysEvaluate(xp, yc);
478 loc[2] = xmap->PhysEvaluate(xp, zc);
480 Array<OneD, NekDouble> uv(2);
481 s->ProjectTo(loc,uv);
487 vector<NodeSharedPtr> honodes;
488 for(
int j = 0; j < niq; j++)
490 Array<OneD, NekDouble> loc;
493 Node(0,loc[0],loc[1],loc[2]));
494 nn->SetCADSurf(surf, s, uvi[j]);
495 honodes.push_back(nn);
498 f->m_faceNodes = honodes;
#define ASSERTL0(condition, msg)
boost::shared_ptr< OptiEdge > OptiEdgeSharedPtr
int PrintProgressbar(const int position, const int goal, const string message, int lastprogress=-1)
Prints a progressbar.
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
MeshSharedPtr m_mesh
Mesh object.
1D Evenly-spaced points using Lagrange polynomial
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
NekMatrix< NekDouble, StandardMatrixTag > DNekMat
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< CADObject > CADObjectSharedPtr
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
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
boost::shared_ptr< Geometry > GeometrySharedPtr
1D Gauss-Lobatto-Legendre quadrature points
2D Nodal Electrostatic Points on a Triangle
boost::shared_ptr< CADCurve > CADCurveSharedPtr