52 namespace NekMeshUtils
57 HOSurfaceMesh::create,
58 "Generates a high-order surface mesh based on CAD");
63 ConfigOption(
true,
"0",
"Perform edge node optimisation.");
73 cout << endl <<
"\tHigh-Order Surface meshing" << endl;
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);
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;
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++)
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++)
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++)
269 nn->SetCADSurf(s[m].first->GetId(), s[m].first, uv);
279 uvb = e->m_n1->GetCADSurfInfo(surf);
280 uve = e->m_n2->GetCADSurfInfo(surf);
283 for (
int k = 0; k < nq; k++)
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;
297 for (
int k = 0; k < nq; k++)
299 all[k * 2 + 0] = uvi[k][0];
300 all[k * 2 + 1] = uvi[k][1];
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++)
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++)
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;
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++)
439 loc[0] = xmap->PhysEvaluate(xp, xc);
440 loc[1] = xmap->PhysEvaluate(xp, yc);
441 loc[2] = xmap->PhysEvaluate(xp, zc);
444 s->ProjectTo(loc,uv);
449 vector<NodeSharedPtr> honodes;
450 for(
int j = 0; j < ni; j++)
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++)
476 loc[0] = xmap->PhysEvaluate(xp, xc);
477 loc[1] = xmap->PhysEvaluate(xp, yc);
478 loc[2] = xmap->PhysEvaluate(xp, zc);
481 s->ProjectTo(loc,uv);
487 vector<NodeSharedPtr> honodes;
488 for(
int j = 0; j < niq; j++)
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
base class for CAD curves.
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.
int GetId()
Return ID of the vertex.
Represents a face comprised of three or more edges.
pair< ModuleType, string > ModuleKey
MeshSharedPtr m_mesh
Mesh object.
1D Evenly-spaced points using Lagrange polynomial
Represents a command-line configuration option.
base class for a cad surface
boost::shared_ptr< Node > NodeSharedPtr
PointsManagerT & PointsManager(void)
Defines a specification for a set of points.
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
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Abstract base class for processing modules.
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
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.