42 namespace NekMeshUtils
49 m_bounds = m_cadcurve->GetBounds();
50 m_curvelength = m_cadcurve->GetTotLength();
52 int(m_curvelength / m_mesh->m_octree->GetMinDelta()) + 10;
53 ds = m_curvelength / (m_numSamplePoints - 1);
59 for (
int i = 0; i < m_numSamplePoints - 1; i++)
61 Ae += ds * (1.0 / m_dst[i][0] + 1.0 / m_dst[i + 1][0]) / 2.0;
70 meshsvalue[1] = m_curvelength;
78 meshsvalue.resize(Ne + 1);
80 meshsvalue[Ne] = m_curvelength;
82 for (
int i = 1; i < Ne; i++)
84 int iterationcounter = 0;
92 NekDouble rhs = EvaluateDS(ski) / Ae * (EvaluatePS(ski) - k);
95 if (abs(lastSki - ski) < 1E-8)
100 ASSERTL0(iterationcounter < 1000000,
"iteration failed");
110 vector<CADVertSharedPtr> verts = m_cadcurve->GetVertex();
111 vector<pair<CADSurfSharedPtr, CADOrientation::Orientation> > s = m_cadcurve->GetAdjSurf();
115 n->SetCADCurve(m_id, m_cadcurve, t);
117 for (
int j = 0; j < s.size(); j++)
119 if (verts[0]->IsDegen() == s[j].first->GetId())
127 n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
129 m_meshpoints.push_back(n);
131 for (
int i = 1; i < meshsvalue.size() - 1; i++)
133 t = m_cadcurve->tAtArcLength(meshsvalue[i]);
134 loc = m_cadcurve->P(t);
136 new Node(m_mesh->m_numNodes++, loc[0], loc[1], loc[2]));
137 n2->SetCADCurve(m_id, m_cadcurve, t);
138 for (
int j = 0; j < s.size(); j++)
141 n2->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
143 m_meshpoints.push_back(n2);
146 n = verts[1]->GetNode();
148 n->SetCADCurve(m_id, m_cadcurve, t);
150 for (
int j = 0; j < s.size(); j++)
152 if (verts[1]->IsDegen() == s[j].first->GetId())
160 n->SetCADSurf(s[j].first->GetId(), s[j].first, uv);
162 m_meshpoints.push_back(n);
164 ASSERTL0(Ne + 1 == m_meshpoints.size(),
165 "incorrect number of points in curve mesh");
168 for (
int i = 0; i < m_meshpoints.size() - 1; i++)
171 new Edge(m_meshpoints[i], m_meshpoints[i + 1]));
172 e->m_parentCAD = m_cadcurve;
173 m_mesh->m_edgeSet.insert(e);
174 m_meshedges.push_back(e);
177 if (m_mesh->m_verbose)
182 cout << scientific <<
"\r\t\tCurve " << m_id << endl
183 <<
"\t\t\tLength: " << m_curvelength << endl
184 <<
"\t\t\tNodes: " << m_meshpoints.size() << endl
185 <<
"\t\t\tSample points: " << m_numSamplePoints << endl
190 void CurveMesh::GetPhiFunction()
192 m_ps.resize(m_numSamplePoints);
193 vector<NekDouble> newPhi;
203 for (
int i = 1; i < m_numSamplePoints; i++)
205 runningInt += (1.0 / m_dst[i - 1][0] + 1.0 / m_dst[i][0]) / 2.0 * ds;
206 newPhi[0] = Ne / Ae * runningInt;
207 newPhi[1] = m_dst[i][1];
217 ASSERTL1(!(s < 0) && !(s > m_curvelength),
"s out of bounds");
223 else if (s == m_curvelength)
225 return m_dst[m_numSamplePoints - 1][0];
228 for (
int i = 0; i < m_numSamplePoints - 1; i++)
230 if (m_dst[i][1] < s && m_dst[i + 1][1] >= s)
246 ASSERTL0(m * s + c == m * s + c,
"DS");
256 ASSERTL1(!(s < 0) && !(s > m_curvelength),
"s out of bounds");
262 else if (s == m_curvelength)
264 return m_ps[m_numSamplePoints - 1][0];
267 for (
int i = 0; i < m_numSamplePoints - 1; i++)
269 if (m_ps[i][1] < s && m_ps[i + 1][1] >= s)
280 cout << a <<
" " << b << endl;
293 ASSERTL0(m * s + c == m * s + c,
"PS");
298 void CurveMesh::GetSampleFunction()
300 m_dst.resize(m_numSamplePoints);
302 vector<NekDouble> dsti;
305 for (
int i = 0; i < m_numSamplePoints; i++)
308 NekDouble t = m_cadcurve->tAtArcLength(dsti[1]);
348 dsti[0] = m_mesh->m_octree->Query(loc);
#define ASSERTL0(condition, msg)
Represents an edge which joins two points.
boost::shared_ptr< Node > NodeSharedPtr
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...