46 namespace NekMeshUtils
61 for (
int i = 0; i < orderedLoops.size(); i++)
63 numPoints += orderedLoops[i].size();
67 ss <<
"3 points required for triangulation, " << numPoints <<
" in loop"
70 for (
int i = 0; i < m_edgeloops.size(); i++)
72 for (
int j = 0; j < m_edgeloops[i].edges.size(); j++)
74 ss << m_edgeloops[i].edges[j]->GetId() <<
" ";
84 vector<Array<OneD, NekDouble> > centers;
85 for (
int i = 0; i < m_edgeloops.size(); i++)
87 centers.push_back(m_edgeloops[i].center);
90 pplanemesh->Assign(orderedLoops, centers, m_id, m_str);
94 pplanemesh->Extract(m_connec);
107 pplanemesh->AssignStiener(m_stienerpoints);
109 pplanemesh->Extract(m_connec);
119 for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
121 (*eit)->m_elLink.clear();
126 for (
int i = 0; i < m_localElements.size(); i++)
128 vector<EdgeSharedPtr> e = m_localElements[i]->GetEdgeList();
129 for (
int j = 0; j < e.size(); j++)
131 e[j]->m_elLink.clear();
133 m_mesh->m_element[m_mesh->m_expDim].push_back(m_localElements[i]);
138 cout << scientific <<
"\r\t\tFace " << m_id << endl
139 <<
"\t\t\tNodes: " << m_localNodes.size() << endl
140 <<
"\t\t\tEdges: " << m_localEdges.size() << endl
141 <<
"\t\t\tTriangles: " << m_localElements.size() << endl
142 <<
"\t\t\tLoops: " << m_edgeloops.size() << endl
143 <<
"\t\t\tSTR: " << m_str << endl
147 void FaceMesh::MakeBL()
149 for (
int i = 1; i < orderedLoops.size(); i++)
151 for (
int j = 0; j < orderedLoops[i].size(); j++)
155 vector<pair<int, CADSurfSharedPtr> > surfs =
156 orderedLoops[i][j]->GetCADSurfs();
158 "point does not have enough surfs to make quad");
162 for (
int s = 0; s < surfs.size(); s++)
164 if (surfs[s].first == m_id)
168 orderedLoops[i][j]->GetCADSurfInfo(surfs[s].first);
170 for (
int k = 0; k < 3; k++)
175 NekDouble mag = sqrt(AN[0] * AN[0] + AN[1] * AN[1] + AN[2] * AN[2]);
176 for (
int k = 0; k < 3; k++)
181 for (
int k = 0; k < 3; k++)
182 tp[k] = m_bl * AN[k] + loc[k];
185 m_cadsurf->ProjectTo(tp, uv);
187 new Node(m_mesh->m_numNodes++, tp[0], tp[1], tp[2]));
188 nn->SetCADSurf(m_id, m_cadsurf, uv);
191 pair<NodeSharedPtr, NodeSharedPtr>(orderedLoops[i][j], nn));
194 orderedLoops[i][j] = nn;
199 void FaceMesh::OptimiseLocalMesh()
210 void FaceMesh::Smoothing()
215 map<int, vector<EdgeSharedPtr> > connectingedges;
217 map<int, vector<ElementSharedPtr> > connectingelements;
219 for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
221 connectingedges[(*eit)->m_n1->m_id].push_back(*eit);
222 connectingedges[(*eit)->m_n2->m_id].push_back(*eit);
225 for (
int i = 0; i < m_localElements.size(); i++)
227 vector<NodeSharedPtr> v = m_localElements[i]->GetVertexList();
228 for (
int j = 0; j < 3; j++)
230 connectingelements[v[j]->m_id].push_back(m_localElements[i]);
235 for (
int q = 0; q < 4; q++)
237 for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
239 if ((*nit)->GetNumCadCurve() > 0)
242 vector<NodeSharedPtr> connodes;
245 vector<EdgeSharedPtr> edges = connectingedges[(*nit)->m_id];
246 vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];
249 for (
int i = 0; i < els.size(); i++)
261 vector<NodeSharedPtr> nodesystem;
262 vector<NekDouble> lamp;
264 for (
int i = 0; i < edges.size(); i++)
266 vector<NekDouble> lambda;
269 if (*nit == edges[i]->m_n1)
271 else if (*nit == edges[i]->m_n2)
274 ASSERTL0(
false,
"could not find node");
279 for (
int j = 0; j < els.size(); j++)
281 vector<NodeSharedPtr> v = els[j]->GetVertexList();
282 if (v[0] == J || v[1] == J || v[2] == J)
289 vector<EdgeSharedPtr> es = els[j]->GetEdgeList();
290 for (
int k = 0; k < es.size(); k++)
292 if (!(es[k]->m_n1 == *nit || es[k]->m_n2 == *nit))
299 ASSERTL0(found,
"failed to find edge to test");
304 NekDouble lam = ((A[0] - uj[0]) * (B[1] - A[1]) -
305 (A[1] - uj[1]) * (B[0] - A[0])) /
306 ((ui[0] - uj[0]) * (B[1] - A[1]) -
307 (ui[1] - uj[1]) * (B[0] - A[0]));
309 if (!(lam < 0) && !(lam > 1))
310 lambda.push_back(lam);
313 if (lambda.size() > 0)
315 sort(lambda.begin(), lambda.end());
318 ud[0] = uj[0] + lambda[0] * (ui[0] - uj[0]);
319 ud[1] = uj[1] + lambda[0] * (ui[1] - uj[1]);
322 new Node(0, locd[0], locd[1], locd[2]));
323 dn->SetCADSurf(m_id, m_cadsurf, ud);
325 nodesystem.push_back(dn);
326 lamp.push_back(lambda[0]);
330 nodesystem.push_back(J);
339 for (
int i = 0; i < nodesystem.size(); i++)
342 ui[0] += uj[0] / nodesystem.size();
343 ui[1] += uj[1] / nodesystem.size();
353 if (!(uvn[0] < bounds[0] || uvn[0] > bounds[1] ||
354 uvn[1] < bounds[2] || uvn[1] > bounds[3]))
357 (*nit)->Move(l2, m_id, uvn);
363 void FaceMesh::DiagonalSwap()
367 map<int, int> idealConnec;
368 map<int, int> actualConnec;
369 map<int, vector<EdgeSharedPtr> > nodetoedge;
372 for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
374 nodetoedge[(*eit)->m_n1->m_id].push_back(*eit);
375 nodetoedge[(*eit)->m_n2->m_id].push_back(*eit);
378 for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
380 if ((*nit)->GetNumCadCurve() == 0)
383 idealConnec[(*nit)->m_id] = 6;
389 vector<NodeSharedPtr> ns;
390 vector<EdgeSharedPtr> e = nodetoedge[(*nit)->m_id];
391 for (
int i = 0; i < e.size(); i++)
397 if (e[i]->m_n1 == (*nit))
398 ns.push_back(e[i]->m_n2);
400 ns.push_back(e[i]->m_n1);
403 "failed to find 2 nodes in the angle system");
405 idealConnec[(*nit)->m_id] =
406 ceil((*nit)->Angle(ns[0], ns[1]) / 3.142 * 3) + 1;
409 for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
411 actualConnec[(*nit)->m_id] = nodetoedge[(*nit)->m_id].size();
416 for (
int q = 0; q < 4; q++)
418 int edgesStart = m_localEdges.size();
420 m_localEdges.clear();
422 int swappedEdges = 0;
426 for (it = edges.begin(); it != edges.end(); it++)
430 if (e->m_elLink.size() != 2)
432 m_localEdges.insert(e);
435 if (e->m_elLink[0].first->GetConf().m_e ==
437 e->m_elLink[1].first->GetConf().m_e ==
440 m_localEdges.insert(e);
450 vector<NodeSharedPtr> nt = tri1->GetVertexList();
454 if (nt[0] != n1 && nt[0] != n2)
460 else if (nt[1] != n1 && nt[1] != n2)
466 else if (nt[2] != n1 && nt[2] != n2)
474 ASSERTL0(
false,
"failed to identify verticies in tri1");
477 nt = tri2->GetVertexList();
479 if (nt[0] != n1 && nt[0] != n2)
483 else if (nt[1] != n1 && nt[1] != n2)
487 else if (nt[2] != n1 && nt[2] != n2)
493 ASSERTL0(
false,
"failed to identify verticies in tri2");
498 ai = A->GetCADSurfInfo(m_id);
499 bi = B->GetCADSurfInfo(m_id);
500 ci = C->GetCADSurfInfo(m_id);
501 di = D->GetCADSurfInfo(m_id);
505 CDA = 0.5 * (-di[0] * ci[1] + ai[0] * ci[1] + ci[0] * di[1] -
506 ai[0] * di[1] - ci[0] * ai[1] + di[0] * ai[1]);
508 CBD = 0.5 * (-bi[0] * ci[1] + di[0] * ci[1] + ci[0] * bi[1] -
509 di[0] * bi[1] - ci[0] * di[1] + bi[0] * di[1]);
513 if (!(CDA > 0.001 && CBD > 0.001))
515 m_localEdges.insert(e);
523 int nodedefectbefore = 0;
525 abs(actualConnec[A->m_id] - idealConnec[A->m_id]);
527 abs(actualConnec[B->m_id] - idealConnec[B->m_id]);
529 abs(actualConnec[C->m_id] - idealConnec[C->m_id]);
531 abs(actualConnec[D->m_id] - idealConnec[D->m_id]);
533 int nodedefectafter = 0;
535 abs(actualConnec[A->m_id] - 1 - idealConnec[A->m_id]);
537 abs(actualConnec[B->m_id] - 1 - idealConnec[B->m_id]);
539 abs(actualConnec[C->m_id] + 1 - idealConnec[C->m_id]);
541 abs(actualConnec[D->m_id] + 1 - idealConnec[D->m_id]);
543 if (nodedefectafter < nodedefectbefore)
550 NekDouble minanglebefore = C->Angle(A, B);
551 minanglebefore = min(minanglebefore, A->Angle(B, C));
552 minanglebefore = min(minanglebefore, B->Angle(A, C));
553 minanglebefore = min(minanglebefore, B->Angle(A, D));
554 minanglebefore = min(minanglebefore, A->Angle(B, D));
555 minanglebefore = min(minanglebefore, D->Angle(A, B));
557 NekDouble minangleafter = C->Angle(B, D);
558 minangleafter = min(minangleafter, D->Angle(B, C));
559 minangleafter = min(minangleafter, B->Angle(C, D));
560 minangleafter = min(minangleafter, C->Angle(A, D));
561 minangleafter = min(minangleafter, A->Angle(C, D));
562 minangleafter = min(minangleafter, D->Angle(A, C));
564 if (minangleafter > minanglebefore)
572 actualConnec[A->m_id]--;
573 actualConnec[B->m_id]--;
574 actualConnec[C->m_id]++;
575 actualConnec[D->m_id]++;
579 CAt = boost::shared_ptr<Edge>(
new Edge(C, A));
580 ADt = boost::shared_ptr<Edge>(
new Edge(A, D));
581 DBt = boost::shared_ptr<Edge>(
new Edge(D, B));
582 BCt = boost::shared_ptr<Edge>(
new Edge(B, C));
584 vector<EdgeSharedPtr> es = tri1->GetEdgeList();
585 for (
int i = 0; i < 3; i++)
596 es = tri2->GetEdgeList();
597 for (
int i = 0; i < 3; i++)
610 vector<pair<ElementSharedPtr, int> > links;
612 links = CA->m_elLink;
613 CA->m_elLink.clear();
614 for (
int i = 0; i < links.size(); i++)
616 if (links[i].first->GetId() == tri1->GetId())
618 CA->m_elLink.push_back(links[i]);
621 links = BC->m_elLink;
622 BC->m_elLink.clear();
623 for (
int i = 0; i < links.size(); i++)
625 if (links[i].first->GetId() == tri1->GetId())
627 BC->m_elLink.push_back(links[i]);
630 links = AD->m_elLink;
631 AD->m_elLink.clear();
632 for (
int i = 0; i < links.size(); i++)
634 if (links[i].first->GetId() == tri2->GetId())
636 AD->m_elLink.push_back(links[i]);
639 links = DB->m_elLink;
640 DB->m_elLink.clear();
641 for (
int i = 0; i < links.size(); i++)
643 if (links[i].first->GetId() == tri2->GetId())
645 DB->m_elLink.push_back(links[i]);
650 vector<NodeSharedPtr> t1, t2;
659 vector<int> tags = tri1->GetTagList();
661 int id1 = tri1->GetId();
662 int id2 = tri2->GetId();
666 tags = tri2->GetTagList();
672 ntri1->CADSurfId = m_id;
673 ntri2->CADSurfId = m_id;
675 vector<EdgeSharedPtr> t1es = ntri1->GetEdgeList();
676 for (
int i = 0; i < 3; i++)
680 ntri1->SetEdge(i, DB);
681 DB->m_elLink.push_back(
682 pair<ElementSharedPtr, int>(ntri1, i));
684 else if (t1es[i] == BC)
686 ntri1->SetEdge(i, BC);
687 BC->m_elLink.push_back(
688 pair<ElementSharedPtr, int>(ntri1, i));
690 else if (t1es[i] == newe)
692 ntri1->SetEdge(i, newe);
693 newe->m_elLink.push_back(
694 pair<ElementSharedPtr, int>(ntri1, i));
698 ASSERTL0(
false,
"weird edge in new tri 1");
701 vector<EdgeSharedPtr> t2es = ntri2->GetEdgeList();
702 for (
int i = 0; i < 3; i++)
706 ntri2->SetEdge(i, CA);
707 CA->m_elLink.push_back(
708 pair<ElementSharedPtr, int>(ntri2, i));
710 else if (t2es[i] == AD)
712 ntri2->SetEdge(i, AD);
713 AD->m_elLink.push_back(
714 pair<ElementSharedPtr, int>(ntri2, i));
716 else if (t2es[i] == newe)
718 ntri2->SetEdge(i, newe);
719 newe->m_elLink.push_back(
720 pair<ElementSharedPtr, int>(ntri2, i));
724 ASSERTL0(
false,
"weird edge in new tri 2");
728 m_localEdges.insert(newe);
730 m_localElements[id1] = ntri1;
731 m_localElements[id2] = ntri2;
737 m_localEdges.insert(e);
741 ASSERTL0(m_localEdges.size() == edgesStart,
"mismatch edge count");
745 void FaceMesh::BuildLocalMesh()
751 int tricomp = m_mesh->m_numcomp++;
756 int quadcomp = m_mesh->m_numcomp++;
757 for (
int i = 0; i < blpairs.size() - 1;
761 vector<NodeSharedPtr> ns;
762 ns.push_back(blpairs[i].first);
763 ns.push_back(blpairs[i].second);
764 ns.push_back(blpairs[i + 1].second);
765 ns.push_back(blpairs[i + 1].first);
767 tags.push_back(quadcomp);
771 m_localElements.push_back(E);
775 vector<NodeSharedPtr> ns;
776 ns.push_back(blpairs.back().first);
777 ns.push_back(blpairs.back().second);
778 ns.push_back(blpairs[0].second);
779 ns.push_back(blpairs[0].first);
781 tags.push_back(quadcomp);
785 m_localElements.push_back(E);
787 for (
int i = 0; i < m_localElements.size(); i++)
789 vector<NodeSharedPtr> nods = m_localElements[i]->GetVertexList();
790 for (
int j = 0; j < nods.size(); j++)
793 m_localNodes.insert(nods[j]);
795 vector<EdgeSharedPtr> edgs = m_localElements[i]->GetEdgeList();
796 for (
int j = 0; j < edgs.size(); j++)
800 if (!(s == m_mesh->m_edgeSet.end()))
803 m_localElements[i]->SetEdge(j, edgs[j]);
806 pair<EdgeSet::iterator, bool> test =
807 m_localEdges.insert(edgs[j]);
812 ->m_elLink.push_back(
813 pair<ElementSharedPtr, int>(m_localElements[i], j));
817 m_localElements[i]->SetEdge(j, *test.first);
819 ->m_elLink.push_back(
820 pair<ElementSharedPtr, int>(m_localElements[i], j));
823 m_localElements[i]->SetId(i);
827 for (
int i = 0; i < m_connec.size(); i++)
832 tags.push_back(tricomp);
837 vector<NodeSharedPtr> nods = E->GetVertexList();
838 for (
int j = 0; j < nods.size(); j++)
841 m_localNodes.insert(nods[j]);
843 vector<EdgeSharedPtr> edgs = E->GetEdgeList();
844 for (
int j = 0; j < edgs.size(); j++)
848 if (!(s == m_mesh->m_edgeSet.end()))
851 E->SetEdge(j, edgs[j]);
854 pair<EdgeSet::iterator, bool> test = m_localEdges.insert(edgs[j]);
859 ->m_elLink.push_back(pair<ElementSharedPtr, int>(E, j));
863 E->SetEdge(j, *test.first);
865 ->m_elLink.push_back(pair<ElementSharedPtr, int>(E, j));
868 E->SetId(m_localElements.size());
869 m_localElements.push_back(E);
873 void FaceMesh::Stretching()
879 NekDouble dxu = int(bnds[1] - bnds[0] < bnds[3] - bnds[2]
881 : (bnds[1] - bnds[0]) / (bnds[3] - bnds[2]) * 40);
882 NekDouble dxv = int(bnds[3] - bnds[2] < bnds[1] - bnds[0]
884 : (bnds[3] - bnds[2]) / (bnds[1] - bnds[0]) * 40);
886 NekDouble du = (bnds[1] - bnds[0]) / dxu;
887 NekDouble dv = (bnds[3] - bnds[2]) / dxv;
891 for (
int i = 0; i < dxu; i++)
893 for (
int j = 0; j < dxv; j++)
896 uv[0] = bnds[0] + i * du;
897 uv[1] = bnds[2] + j * dv;
904 NekDouble ru = sqrt(r[3] * r[3] + r[4] * r[4] + r[5] * r[5]);
905 NekDouble rv = sqrt(r[6] * r[6] + r[7] * r[7] + r[8] * r[8]);
921 bool FaceMesh::Validate()
926 int pointBefore = m_stienerpoints.size();
927 for (
int i = 0; i < m_connec.size(); i++)
933 r[0] = m_connec[i][0]->Distance(m_connec[i][1]);
934 r[1] = m_connec[i][1]->Distance(m_connec[i][2]);
935 r[2] = m_connec[i][2]->Distance(m_connec[i][0]);
937 triDelta[0] = m_octree->Query(m_connec[i][0]->GetLoc());
938 triDelta[1] = m_octree->Query(m_connec[i][1]->GetLoc());
939 triDelta[2] = m_octree->Query(m_connec[i][2]->GetLoc());
943 if (r[0] < triDelta[0] && r[2] < triDelta[0])
946 if (r[1] < triDelta[1] && r[0] < triDelta[1])
949 if (r[2] < triDelta[2] && r[1] < triDelta[2])
955 ainfo = m_connec[i][0]->GetCADSurfInfo(m_id);
956 binfo = m_connec[i][1]->GetCADSurfInfo(m_id);
957 cinfo = m_connec[i][2]->GetCADSurfInfo(m_id);
960 uvc[0] = (ainfo[0] + binfo[0] + cinfo[0]) / 3.0;
961 uvc[1] = (ainfo[1] + binfo[1] + cinfo[1]) / 3.0;
966 if (m_stienerpoints.size() == pointBefore)
983 new Node(m_mesh->m_numNodes++, np[0], np[1], np[2]));
987 for (
int i = 0; i < orderedLoops.size(); i++)
989 for (
int j = 0; j < orderedLoops[i].size(); j++)
991 NekDouble r = orderedLoops[i][j]->Distance(n);
993 if (r < npDelta / 2.0)
1003 for (
int i = 0; i < m_stienerpoints.size(); i++)
1005 NekDouble r = m_stienerpoints[i]->Distance(n);
1007 if (r < npDelta / 2.0)
1017 n->SetCADSurf(m_id, m_cadsurf, uv);
1018 m_stienerpoints.push_back(n);
1022 void FaceMesh::OrientateCurves()
1025 for (
int i = 0; i < m_edgeloops.size(); i++)
1027 vector<NodeSharedPtr> cE;
1028 for (
int j = 0; j < m_edgeloops[i].edges.size(); j++)
1030 int cid = m_edgeloops[i].edges[j]->GetId();
1031 vector<NodeSharedPtr> edgePoints =
1032 m_curvemeshes[cid]->GetMeshPoints();
1034 int numPoints = m_curvemeshes[cid]->GetNumPoints();
1036 if (m_edgeloops[i].edgeo[j] == 0)
1038 for (
int k = 0; k < numPoints - 1; k++)
1040 cE.push_back(edgePoints[k]);
1045 for (
int k = numPoints - 1; k > 0; k--)
1047 cE.push_back(edgePoints[k]);
1051 orderedLoops.push_back(cE);
1055 for (
int i = 0; i < orderedLoops.size(); i++)
1058 for (
int j = 0; j < orderedLoops[i].size() - 1; j++)
1061 n1info = orderedLoops[i][j]->GetCADSurfInfo(m_id);
1062 n2info = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
1064 area += -n2info[1] * (n2info[0] - n1info[0]) +
1065 n1info[0] * (n2info[1] - n1info[1]);
1068 m_edgeloops[i].area = area;
1076 for (
int i = 0; i < m_edgeloops.size() - 1; i++)
1078 if (fabs(m_edgeloops[i].area) < fabs(m_edgeloops[i + 1].area))
1081 vector<NodeSharedPtr> orderedlooptmp = orderedLoops[i];
1082 EdgeLoop edgeLoopstmp = m_edgeloops[i];
1084 orderedLoops[i] = orderedLoops[i + 1];
1085 m_edgeloops[i] = m_edgeloops[i + 1];
1087 orderedLoops[i + 1] = orderedlooptmp;
1088 m_edgeloops[i + 1] = edgeLoopstmp;
1096 for (
int i = 0; i < orderedLoops.size(); i++)
1100 n1 = orderedLoops[i][0];
1101 n2 = orderedLoops[i][1];
1104 n1info = n1->GetCADSurfInfo(m_id);
1105 n2info = n2->GetCADSurfInfo(m_id);
1108 NekDouble mag = sqrt((n1info[0] - n2info[0]) * (n1info[0] - n2info[0]) +
1109 (n1info[1] - n2info[1]) * (n1info[1] - n2info[1]));
1111 N[0] = -1.0 * (n2info[1] - n1info[1]) / mag;
1112 N[1] = (n2info[0] - n1info[0]) / mag;
1115 P[0] = (n1info[0] + n2info[0]) / 2.0 + 1e-6 * N[0];
1116 P[1] = (n1info[1] + n2info[1]) / 2.0 + 1e-6 * N[1];
1121 for (
int j = 0; j < orderedLoops[i].size() - 1; j++)
1124 nt1 = orderedLoops[i][j]->GetCADSurfInfo(m_id);
1125 nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
1127 if (fabs(nt2[1] - nt1[1]) < 1e-30)
1130 NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1131 NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1133 if (!(lam < 0) && !(lam > 1) && S > 0)
1140 nt1 = orderedLoops[i].back()->GetCADSurfInfo(m_id);
1141 nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);
1143 if (fabs(nt2[1] - nt1[1]) < 1e-30)
1146 NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1147 NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1149 if (!(lam < 0) && !(lam > 1) && S > 0)
1154 if (intercepts % 2 == 0)
1156 P[0] = (n1info[0] + n2info[0]) / 2.0 - 1e-6 * N[0];
1157 P[1] = (n1info[1] + n2info[1]) / 2.0 - 1e-6 * N[1];
1159 for (
int j = 0; j < orderedLoops[i].size() - 1; j++)
1162 nt1 = orderedLoops[i][j]->GetCADSurfInfo(m_id);
1163 nt2 = orderedLoops[i][j + 1]->GetCADSurfInfo(m_id);
1165 if (fabs(nt2[1] - nt1[1]) < 1e-30)
1168 NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1169 NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1171 if (!(lam < 0) && !(lam > 1) && S > 0)
1178 nt1 = orderedLoops[i].back()->GetCADSurfInfo(m_id);
1179 nt2 = orderedLoops[i][0]->GetCADSurfInfo(m_id);
1181 if (fabs(nt2[1] - nt1[1]) < 1e-30)
1184 NekDouble lam = (P[1] - nt1[1]) / (nt2[1] - nt1[1]);
1185 NekDouble S = nt1[0] - P[0] + (nt2[0] - nt1[0]) * lam;
1187 if (!(lam < 0) && !(lam > 1) && S > 0)
1192 if (intercepts % 2 == 0)
1194 cerr <<
"still failed to find point inside loop" << endl;
1198 m_edgeloops[i].center = P;
1201 if (m_edgeloops[0].area < 0)
1203 vector<NodeSharedPtr> tmp = orderedLoops[0];
1204 reverse(tmp.begin(), tmp.end());
1205 orderedLoops[0] = tmp;
1208 for (
int i = 1; i < orderedLoops.size(); i++)
1210 if (m_edgeloops[i].area > 0)
1212 vector<NodeSharedPtr> tmp = orderedLoops[i];
1213 reverse(tmp.begin(), tmp.end());
1214 orderedLoops[i] = tmp;
#define ASSERTL0(condition, msg)
Basic information about an element.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Represents an edge which joins two points.
boost::shared_ptr< TriangleInterface > TriangleInterfaceSharedPtr
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
struct which descibes a collection of cad edges which for a loop on the cad surface ...
ElementFactory & GetElementFactory()
Represents a point in the domain.
boost::shared_ptr< Node > NodeSharedPtr
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Element > ElementSharedPtr
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet