44 namespace NekMeshUtils
47 bool FaceMesh::ValidateCurves()
49 vector<int> curvesInSurface;
50 for (
int i = 0; i < m_edgeloops.size(); i++)
52 for (
int j = 0; j < m_edgeloops[i]->edges.size(); j++)
54 curvesInSurface.push_back(m_edgeloops[i]->edges[j]->GetId());
60 for (
int i = 0; i < curvesInSurface.size(); i++)
62 vector<EdgeSharedPtr> es =
63 m_curvemeshes[curvesInSurface[i]]->GetMeshEdges();
65 for (
int j = i; j < curvesInSurface.size(); j++)
72 vector<EdgeSharedPtr> es2 =
73 m_curvemeshes[curvesInSurface[j]]->GetMeshEdges();
75 for (
int l = 0; l < es.size(); l++)
79 for (
int k = 0; k < es2.size(); k++)
81 if (es[l]->m_n1 == es2[k]->m_n1 ||
82 es[l]->m_n1 == es2[k]->m_n2 ||
83 es[l]->m_n2 == es2[k]->m_n1 ||
84 es[l]->m_n2 == es2[k]->m_n2)
90 es2[k]->m_n1->GetCADSurfInfo(m_id);
92 es2[k]->m_n2->GetCADSurfInfo(m_id);
94 NekDouble den = (P4[0] - P3[0]) * (P2[1] - P1[1]) -
95 (P2[0] - P1[0]) * (P4[1] - P3[1]);
102 NekDouble t = ((P1[0] - P3[0]) * (P4[1] - P3[1]) -
103 (P4[0] - P3[0]) * (P1[1] - P3[1])) /
106 (P1[0] - P3[0] + t * (P2[0] - P1[0])) / (P4[0] - P3[0]);
108 if (t < 1.0 && t > 0.0 && u < 1.0 && u > 0.0)
111 uv[0] = P1[0] + t * (P2[0] - P1[0]);
112 uv[1] = P1[1] + t * (P2[1] - P1[1]);
115 <<
"Curve mesh error at " << loc[0] <<
" "
116 << loc[1] <<
" " << loc[2] <<
" on face " << m_id
133 for (
int i = 0; i < orderedLoops.size(); i++)
135 numPoints += orderedLoops[i].size();
136 for (
int j = 0; j < orderedLoops[i].size(); j++)
138 m_inBoundary.insert(orderedLoops[i][j]);
143 ss <<
"3 points required for triangulation, " << numPoints <<
" in loop"
146 for (
int i = 0; i < m_edgeloops.size(); i++)
148 for (
int j = 0; j < m_edgeloops[i]->edges.size(); j++)
150 ss << m_edgeloops[i]->edges[j]->GetId() <<
" ";
160 vector<Array<OneD, NekDouble> > centers;
161 for (
int i = 0; i < m_edgeloops.size(); i++)
163 centers.push_back(m_edgeloops[i]->center);
166 pplanemesh->Assign(orderedLoops, centers, m_id, m_str);
170 pplanemesh->Extract(m_connec);
185 pplanemesh->AssignStiener(m_stienerpoints);
187 pplanemesh->Extract(m_connec);
200 for (
int i = 0; i < m_localElements.size(); i++)
202 vector<EdgeSharedPtr> e = m_localElements[i]->GetEdgeList();
203 for (
int j = 0; j < e.size(); j++)
205 e[j]->m_elLink.clear();
207 m_mesh->m_element[2].push_back(m_localElements[i]);
210 if (m_mesh->m_verbose)
215 cout << scientific <<
"\r\t\tFace " << m_id << endl
216 <<
"\t\t\tNodes: " << m_localNodes.size() << endl
217 <<
"\t\t\tEdges: " << m_localEdges.size() << endl
218 <<
"\t\t\tTriangles: " << m_localElements.size() << endl
219 <<
"\t\t\tLoops: " << m_edgeloops.size() << endl
220 <<
"\t\t\tSTR: " << m_str << endl
225 void FaceMesh::OptimiseLocalMesh()
237 void FaceMesh::Smoothing()
244 map<int, vector<EdgeSharedPtr> > connectingedges;
246 map<int, vector<ElementSharedPtr> > connectingelements;
248 for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
250 connectingedges[(*eit)->m_n1->m_id].push_back(*eit);
251 connectingedges[(*eit)->m_n2->m_id].push_back(*eit);
254 for (
int i = 0; i < m_localElements.size(); i++)
256 vector<NodeSharedPtr> v = m_localElements[i]->GetVertexList();
257 for (
int j = 0; j < 3; j++)
259 connectingelements[v[j]->m_id].push_back(m_localElements[i]);
264 for (
int q = 0; q < 4; q++)
266 for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
271 if (f != m_inBoundary.end())
277 vector<NodeSharedPtr> connodes;
279 vector<EdgeSharedPtr> edges = connectingedges[(*nit)->m_id];
280 vector<ElementSharedPtr> els = connectingelements[(*nit)->m_id];
282 vector<NodeSharedPtr> nodesystem;
283 vector<NekDouble> lamp;
285 for (
int i = 0; i < edges.size(); i++)
287 vector<NekDouble> lambda;
290 if (*nit == edges[i]->m_n1)
294 else if (*nit == edges[i]->m_n2)
300 ASSERTL0(
false,
"could not find node");
306 for (
int j = 0; j < els.size(); j++)
308 vector<NodeSharedPtr> v = els[j]->GetVertexList();
311 if (v[0] == J || v[1] == J || v[2] == J)
319 vector<EdgeSharedPtr> es = els[j]->GetEdgeList();
320 for (
int k = 0; k < es.size(); k++)
322 if (!(es[k]->m_n1 == *nit || es[k]->m_n2 == *nit))
329 ASSERTL0(found,
"failed to find edge to test");
334 NekDouble lam = ((A[0] - uj[0]) * (B[1] - A[1]) -
335 (A[1] - uj[1]) * (B[0] - A[0])) /
336 ((ui[0] - uj[0]) * (B[1] - A[1]) -
337 (ui[1] - uj[1]) * (B[0] - A[0]));
339 if (!(lam < 0) && !(lam > 1))
340 lambda.push_back(lam);
343 if (lambda.size() > 0)
345 sort(lambda.begin(), lambda.end());
348 ud[0] = uj[0] + lambda[0] * (ui[0] - uj[0]);
349 ud[1] = uj[1] + lambda[0] * (ui[1] - uj[1]);
352 new Node(0, locd[0], locd[1], locd[2]));
353 dn->SetCADSurf(m_id, m_cadsurf, ud);
355 nodesystem.push_back(dn);
356 lamp.push_back(lambda[0]);
360 nodesystem.push_back(J);
369 for (
int i = 0; i < nodesystem.size(); i++)
372 u0[0] += uj[0] / nodesystem.size();
373 u0[1] += uj[1] / nodesystem.size();
434 bool inbounds =
true;
435 if (u0[0] < bounds[0])
439 else if (u0[0] > bounds[1])
443 else if (u0[1] < bounds[2])
447 else if (u0[1] > bounds[3])
483 (*nit)->Move(l2, m_id, u0);
488 void FaceMesh::DiagonalSwap()
490 map<int, int> idealConnec;
491 map<int, int> actualConnec;
492 map<int, vector<EdgeSharedPtr> > nodetoedge;
495 for (eit = m_localEdges.begin(); eit != m_localEdges.end(); eit++)
497 nodetoedge[(*eit)->m_n1->m_id].push_back(*eit);
498 nodetoedge[(*eit)->m_n2->m_id].push_back(*eit);
501 for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
505 if (f == m_inBoundary.end())
508 idealConnec[(*nit)->m_id] = 6;
532 idealConnec[(*nit)->m_id] = 4;
535 for (nit = m_localNodes.begin(); nit != m_localNodes.end(); nit++)
537 actualConnec[(*nit)->m_id] = nodetoedge[(*nit)->m_id].size();
542 for (
int q = 0; q < 4; q++)
544 int edgesStart = m_localEdges.size();
546 m_localEdges.clear();
548 int swappedEdges = 0;
552 for (it = edges.begin(); it != edges.end(); it++)
558 if (f1 != m_inBoundary.end() && f2 != m_inBoundary.end())
560 m_localEdges.insert(e);
570 vector<NodeSharedPtr> nt = tri1->GetVertexList();
574 if (nt[0] != n1 && nt[0] != n2)
580 else if (nt[1] != n1 && nt[1] != n2)
586 else if (nt[2] != n1 && nt[2] != n2)
594 ASSERTL0(
false,
"failed to identify verticies in tri1");
597 nt = tri2->GetVertexList();
599 if (nt[0] != n1 && nt[0] != n2)
603 else if (nt[1] != n1 && nt[1] != n2)
607 else if (nt[2] != n1 && nt[2] != n2)
613 ASSERTL0(
false,
"failed to identify verticies in tri2");
628 ai = A->GetCADSurfInfo(m_id);
629 bi = B->GetCADSurfInfo(m_id);
630 ci = C->GetCADSurfInfo(m_id);
631 di = D->GetCADSurfInfo(m_id);
635 CDA = 0.5 * (-di[0] * ci[1] + ai[0] * ci[1] + ci[0] * di[1] -
636 ai[0] * di[1] - ci[0] * ai[1] + di[0] * ai[1]);
638 CBD = 0.5 * (-bi[0] * ci[1] + di[0] * ci[1] + ci[0] * bi[1] -
639 di[0] * bi[1] - ci[0] * di[1] + bi[0] * di[1]);
643 if (!(CDA > 0.001 && CBD > 0.001))
645 m_localEdges.insert(e);
653 int nodedefectbefore = 0;
655 abs(actualConnec[A->m_id] - idealConnec[A->m_id]);
657 abs(actualConnec[B->m_id] - idealConnec[B->m_id]);
659 abs(actualConnec[C->m_id] - idealConnec[C->m_id]);
661 abs(actualConnec[D->m_id] - idealConnec[D->m_id]);
663 int nodedefectafter = 0;
665 abs(actualConnec[A->m_id] - 1 - idealConnec[A->m_id]);
667 abs(actualConnec[B->m_id] - 1 - idealConnec[B->m_id]);
669 abs(actualConnec[C->m_id] + 1 - idealConnec[C->m_id]);
671 abs(actualConnec[D->m_id] + 1 - idealConnec[D->m_id]);
673 if (nodedefectafter < nodedefectbefore)
680 NekDouble minanglebefore = C->Angle(A, B);
681 minanglebefore = min(minanglebefore, A->Angle(B, C));
682 minanglebefore = min(minanglebefore, B->Angle(A, C));
683 minanglebefore = min(minanglebefore, B->Angle(A, D));
684 minanglebefore = min(minanglebefore, A->Angle(B, D));
685 minanglebefore = min(minanglebefore, D->Angle(A, B));
687 NekDouble minangleafter = C->Angle(B, D);
688 minangleafter = min(minangleafter, D->Angle(B, C));
689 minangleafter = min(minangleafter, B->Angle(C, D));
690 minangleafter = min(minangleafter, C->Angle(A, D));
691 minangleafter = min(minangleafter, A->Angle(C, D));
692 minangleafter = min(minangleafter, D->Angle(A, C));
694 if (minangleafter > minanglebefore)
702 actualConnec[A->m_id]--;
703 actualConnec[B->m_id]--;
704 actualConnec[C->m_id]++;
705 actualConnec[D->m_id]++;
709 CAt = boost::shared_ptr<Edge>(
new Edge(C, A));
710 ADt = boost::shared_ptr<Edge>(
new Edge(A, D));
711 DBt = boost::shared_ptr<Edge>(
new Edge(D, B));
712 BCt = boost::shared_ptr<Edge>(
new Edge(B, C));
714 vector<EdgeSharedPtr> es = tri1->GetEdgeList();
715 for (
int i = 0; i < 3; i++)
726 es = tri2->GetEdgeList();
727 for (
int i = 0; i < 3; i++)
740 vector<pair<ElementSharedPtr, int> > links;
742 links = CA->m_elLink;
743 CA->m_elLink.clear();
744 for (
int i = 0; i < links.size(); i++)
746 if (links[i].first->GetId() == tri1->GetId())
750 CA->m_elLink.push_back(links[i]);
753 links = BC->m_elLink;
754 BC->m_elLink.clear();
755 for (
int i = 0; i < links.size(); i++)
757 if (links[i].first->GetId() == tri1->GetId())
761 BC->m_elLink.push_back(links[i]);
764 links = AD->m_elLink;
765 AD->m_elLink.clear();
766 for (
int i = 0; i < links.size(); i++)
768 if (links[i].first->GetId() == tri2->GetId())
772 AD->m_elLink.push_back(links[i]);
775 links = DB->m_elLink;
776 DB->m_elLink.clear();
777 for (
int i = 0; i < links.size(); i++)
779 if (links[i].first->GetId() == tri2->GetId())
783 DB->m_elLink.push_back(links[i]);
788 vector<NodeSharedPtr> t1, t2;
797 vector<int> tags = tri1->GetTagList();
799 int id1 = tri1->GetId();
800 int id2 = tri2->GetId();
804 tags = tri2->GetTagList();
810 ntri1->m_parentCAD = m_cadsurf;
811 ntri2->m_parentCAD = m_cadsurf;
813 vector<EdgeSharedPtr> t1es = ntri1->GetEdgeList();
814 for (
int i = 0; i < 3; i++)
818 ntri1->SetEdge(i, DB);
819 DB->m_elLink.push_back(
820 pair<ElementSharedPtr, int>(ntri1, i));
822 else if (t1es[i] == BC)
824 ntri1->SetEdge(i, BC);
825 BC->m_elLink.push_back(
826 pair<ElementSharedPtr, int>(ntri1, i));
828 else if (t1es[i] == newe)
830 ntri1->SetEdge(i, newe);
831 newe->m_elLink.push_back(
832 pair<ElementSharedPtr, int>(ntri1, i));
836 ASSERTL0(
false,
"weird edge in new tri 1");
839 vector<EdgeSharedPtr> t2es = ntri2->GetEdgeList();
840 for (
int i = 0; i < 3; i++)
844 ntri2->SetEdge(i, CA);
845 CA->m_elLink.push_back(
846 pair<ElementSharedPtr, int>(ntri2, i));
848 else if (t2es[i] == AD)
850 ntri2->SetEdge(i, AD);
851 AD->m_elLink.push_back(
852 pair<ElementSharedPtr, int>(ntri2, i));
854 else if (t2es[i] == newe)
856 ntri2->SetEdge(i, newe);
857 newe->m_elLink.push_back(
858 pair<ElementSharedPtr, int>(ntri2, i));
862 ASSERTL0(
false,
"weird edge in new tri 2");
866 m_localEdges.insert(newe);
868 m_localElements[id1] = ntri1;
869 m_localElements[id2] = ntri2;
875 m_localEdges.insert(e);
879 ASSERTL0(m_localEdges.size() == edgesStart,
"mismatch edge count");
883 void FaceMesh::BuildLocalMesh()
890 for (
int i = 0; i < m_connec.size(); i++)
895 tags.push_back(m_compId);
898 E->m_parentCAD = m_cadsurf;
900 vector<NodeSharedPtr> nods = E->GetVertexList();
901 for (
int j = 0; j < nods.size(); j++)
904 m_localNodes.insert(nods[j]);
907 E->SetId(m_localElements.size());
908 m_localElements.push_back(E);
911 for (
int i = 0; i < m_localElements.size(); ++i)
913 for (
int j = 0; j < m_localElements[i]->GetEdgeCount(); ++j)
915 pair<EdgeSet::iterator, bool> testIns;
919 if (!(s == m_mesh->m_edgeSet.end()))
922 m_localElements[i]->SetEdge(j, *s);
925 testIns = m_localEdges.insert(ed);
930 ed2->m_elLink.push_back(
931 pair<ElementSharedPtr, int>(m_localElements[i], j));
936 m_localElements[i]->SetEdge(j, e2);
937 e2->m_elLink.push_back(
938 pair<ElementSharedPtr, int>(m_localElements[i], j));
944 void FaceMesh::Stretching()
950 NekDouble dxu = int(bnds[1] - bnds[0] < bnds[3] - bnds[2]
952 : (bnds[1] - bnds[0]) / (bnds[3] - bnds[2]) * 40);
953 NekDouble dxv = int(bnds[3] - bnds[2] < bnds[1] - bnds[0]
955 : (bnds[3] - bnds[2]) / (bnds[1] - bnds[0]) * 40);
957 NekDouble du = (bnds[1] - bnds[0]) / dxu;
958 NekDouble dv = (bnds[3] - bnds[2]) / dxv;
962 for (
int i = 0; i < dxu; i++)
964 for (
int j = 0; j < dxv; j++)
967 uv[0] = bnds[0] + i * du;
968 uv[1] = bnds[2] + j * dv;
979 NekDouble ru = sqrt(r[3] * r[3] + r[4] * r[4] + r[5] * r[5]);
980 NekDouble rv = sqrt(r[6] * r[6] + r[7] * r[7] + r[8] * r[8]);
998 bool FaceMesh::Validate()
1003 int pointBefore = m_stienerpoints.size();
1004 for (
int i = 0; i < m_connec.size(); i++)
1008 vector<Array<OneD, NekDouble> > info;
1010 for (
int j = 0; j < 3; j++)
1012 info.push_back(m_connec[i][j]->GetCADSurfInfo(m_id));
1015 r[0] = m_connec[i][0]->Distance(m_connec[i][1]);
1016 r[1] = m_connec[i][1]->Distance(m_connec[i][2]);
1017 r[2] = m_connec[i][2]->Distance(m_connec[i][0]);
1019 a[0] = m_connec[i][0]->Angle(m_connec[i][1]->GetLoc(),
1020 m_connec[i][2]->GetLoc(),
1021 m_cadsurf->N(info[0]));
1022 a[1] = m_connec[i][1]->Angle(m_connec[i][2]->GetLoc(),
1023 m_connec[i][0]->GetLoc(),
1024 m_cadsurf->N(info[1]));
1025 a[2] = m_connec[i][2]->Angle(m_connec[i][0]->GetLoc(),
1026 m_connec[i][1]->GetLoc(),
1027 m_cadsurf->N(info[2]));
1029 NekDouble d1 = m_mesh->m_octree->Query(m_connec[i][0]->GetLoc());
1030 NekDouble d2 = m_mesh->m_octree->Query(m_connec[i][1]->GetLoc());
1031 NekDouble d3 = m_mesh->m_octree->Query(m_connec[i][2]->GetLoc());
1034 uvc[0] = (info[0][0] + info[1][0] + info[2][0]) / 3.0;
1035 uvc[1] = (info[0][1] + info[1][1] + info[2][1]) / 3.0;
1038 NekDouble d4 = m_mesh->m_octree->Query(locc);
1040 NekDouble d = (d1 + d2 + d3 + d4) / 4.0;
1042 vector<bool> valid(3);
1043 valid[0] = r[0] < d * 1.5;
1044 valid[1] = r[1] < d * 1.5;
1045 valid[2] = r[2] < d * 1.5;
1047 vector<bool> angValid(3);
1048 angValid[0] = a[0] / M_PI * 180.0 > 20.0 && a[0] / M_PI * 180.0 < 120.0;
1049 angValid[1] = a[1] / M_PI * 180.0 > 20.0 && a[1] / M_PI * 180.0 < 120.0;
1050 angValid[2] = a[2] / M_PI * 180.0 > 20.0 && a[2] / M_PI * 180.0 < 120.0;
1053 int numAngValid = 0;
1054 for (
int j = 0; j < 3; j++)
1072 if (numValid != 3 || numAngValid != 3)
1099 if (m_stienerpoints.size() == pointBefore)
1113 NekDouble npDelta = m_mesh->m_octree->Query(np);
1116 new Node(m_mesh->m_numNodes++, np[0], np[1], np[2]));
1120 for (
int i = 0; i < orderedLoops.size(); i++)
1122 for (
int j = 0; j < orderedLoops[i].size(); j++)
1124 NekDouble r = orderedLoops[i][j]->Distance(n);
1126 if (r < npDelta / 2.0)
1136 for (
int i = 0; i < m_stienerpoints.size(); i++)
1138 NekDouble r = m_stienerpoints[i]->Distance(n);
1140 if (r < npDelta / 2.0)
1150 n->SetCADSurf(m_id, m_cadsurf, uv);
1151 m_stienerpoints.push_back(n);
1155 void FaceMesh::OrientateCurves()
1158 orderedLoops.clear();
1161 for (
int i = 0; i < m_edgeloops.size(); i++)
1163 vector<NodeSharedPtr> cE;
1164 for (
int j = 0; j < m_edgeloops[i]->edges.size(); j++)
1166 int cid = m_edgeloops[i]->edges[j]->GetId();
1167 vector<NodeSharedPtr> edgePoints =
1168 m_curvemeshes[cid]->GetMeshPoints();
1170 int numPoints = m_curvemeshes[cid]->GetNumPoints();
1174 for (
int k = 0; k < numPoints - 1; k++)
1176 cE.push_back(edgePoints[k]);
1181 for (
int k = numPoints - 1; k > 0; k--)
1183 cE.push_back(edgePoints[k]);
1187 orderedLoops.push_back(cE);
#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...
ElementFactory & GetElementFactory()
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