41 namespace NekMeshUtils
44 void CurveMesh::Mesh()
46 m_bounds = m_cadcurve->Bounds();
47 m_curvelength = m_cadcurve->GetTotLength();
48 m_numSamplePoints = int(m_curvelength / m_octree->GetMinDelta()) + 5;
49 ds = m_curvelength / (m_numSamplePoints - 1);
55 for (
int i = 0; i < m_numSamplePoints - 1; i++)
57 Ae += ds * (1.0 / m_dst[i][0] + 1.0 / m_dst[i + 1][0]) / 2.0;
66 meshsvalue[1] = m_curvelength;
74 meshsvalue.resize(Ne + 1);
76 meshsvalue[Ne] = m_curvelength;
78 for (
int i = 1; i < Ne; i++)
80 int iterationcounter = 0;
88 NekDouble rhs = EvaluateDS(ski) / Ae * (EvaluatePS(ski) - k);
91 if (abs(lastSki - ski) < 1E-10)
96 ASSERTL0(iterationcounter < 1000000,
"iteration failed");
106 vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
107 vector<CADSurfSharedPtr> s = m_cadcurve->GetAdjSurf();
108 ASSERTL0(s.size() == 2,
"invalid curve");
112 n->SetCADCurve(m_id, m_cadcurve, t);
114 for (
int j = 0; j < 2; j++)
116 if (verts[0]->IsDegen() == s[j]->GetId())
123 n->SetCADSurf(s[j]->GetId(), s[j], uv);
125 m_meshpoints.push_back(n);
127 for (
int i = 1; i < meshsvalue.size() - 1; i++)
129 t = m_cadcurve->tAtArcLength(meshsvalue[i]);
130 loc = m_cadcurve->P(t);
132 new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
133 n2->SetCADCurve(m_id, m_cadcurve, t);
134 for (
int j = 0; j < 2; j++)
137 n2->SetCADSurf(s[j]->GetId(), s[j], uv);
139 m_meshpoints.push_back(n2);
142 n = verts[1]->GetNode();
144 n->SetCADCurve(m_id, m_cadcurve, t);
146 for (
int j = 0; j < 2; j++)
148 if (verts[1]->IsDegen() == s[j]->GetId())
155 n->SetCADSurf(s[j]->GetId(), s[j], uv);
157 m_meshpoints.push_back(n);
159 ASSERTL0(Ne + 1 == m_meshpoints.size(),
160 "incorrect number of points in curve mesh");
162 for (
int i = 0; i < m_meshpoints.size(); i++)
165 for (
int j = 0; j < 2; j++)
168 m_meshpoints[i]->SetCADSurf(s[j]->GetId(), s[j], uv);
220 for (
int i = 0; i < m_meshpoints.size() - 1; i++)
223 new Edge(m_meshpoints[i], m_meshpoints[i + 1]));
224 e->CADCurveId = m_id;
225 e->CADCurve = m_cadcurve;
227 m_mesh->m_edgeSet.insert(e);
232 cout << scientific <<
"\r\t\tCurve " << m_id << endl
233 <<
"\t\t\tLength: " << m_curvelength << endl
234 <<
"\t\t\tNodes: " << m_meshpoints.size() << endl
235 <<
"\t\t\tSample points: " << m_numSamplePoints << endl
239 void CurveMesh::GetPhiFunction()
241 m_ps.resize(m_numSamplePoints);
242 vector<NekDouble> newPhi;
252 for (
int i = 1; i < m_numSamplePoints; i++)
254 runningInt += (1.0 / m_dst[i - 1][0] + 1.0 / m_dst[i][0]) / 2.0 * ds;
255 newPhi[0] = Ne / Ae * runningInt;
256 newPhi[1] = m_dst[i][1];
270 else if (s == m_curvelength)
272 return m_dst[m_numSamplePoints - 1][0];
275 for (
int i = 0; i < m_numSamplePoints - 1; i++)
277 if (m_dst[i][1] < s && m_dst[i + 1][1] >= s)
293 ASSERTL0(m * s + c == m * s + c,
"DS");
307 else if (s == m_curvelength)
309 return m_ps[m_numSamplePoints - 1][0];
312 for (
int i = 0; i < m_numSamplePoints - 1; i++)
314 if (m_ps[i][1] < s && m_ps[i + 1][1] >= s)
325 cout << a <<
" " << b << endl;
338 ASSERTL0(m * s + c == m * s + c,
"PS");
343 void CurveMesh::GetSampleFunction()
345 m_dst.resize(m_numSamplePoints);
348 vector<NekDouble> dsti;
351 dsti[0] = m_octree->Query(loc);
353 dsti[2] = m_bounds[0];
357 for (
int i = 1; i < m_numSamplePoints; i++)
360 NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
362 loc = m_cadcurve->P(t);
364 dsti[0] = m_octree->Query(loc);
#define ASSERTL0(condition, msg)
Represents an edge which joins two points.
Represents a point in the domain.
boost::shared_ptr< Node > NodeSharedPtr
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.