diff --git a/libsrc/general/hashtabl.hpp b/libsrc/general/hashtabl.hpp index 4776c626..66dd934b 100644 --- a/libsrc/general/hashtabl.hpp +++ b/libsrc/general/hashtabl.hpp @@ -717,6 +717,27 @@ public: } }; + + +template +inline ostream & operator<< (ostream & ost, const INDEX_2_CLOSED_HASHTABLE & ht) +{ + for (int i = 0; i < ht.Size(); i++) + if (ht.UsedPos(i)) + { + INDEX_2 hash; + T data; + ht.GetData (i, hash, data); + ost << "hash = " << hash << ", data = " << data << endl; + } + return ost; +} + + + + + + class BASE_INDEX_3_CLOSED_HASHTABLE { protected: @@ -724,7 +745,6 @@ protected: int invalid; protected: - // public: //SZ BASE_INDEX_3_CLOSED_HASHTABLE (int size) : hash(size) { @@ -896,6 +916,20 @@ public: +template +inline ostream & operator<< (ostream & ost, const INDEX_3_CLOSED_HASHTABLE & ht) +{ + for (int i = 0; i < ht.Size(); i++) + if (ht.UsedPos(i)) + { + INDEX_3 hash; + T data; + ht.GetData (i, hash, data); + ost << "hash = " << hash << ", data = " << data << endl; + } + return ost; +} + diff --git a/libsrc/general/profiler.cpp b/libsrc/general/profiler.cpp index f17be951..60efaa08 100644 --- a/libsrc/general/profiler.cpp +++ b/libsrc/general/profiler.cpp @@ -42,19 +42,20 @@ namespace netgen // thus we use the C-variant: - - /* - char filename[100]; + if (getenv ("NGPROFILE")) + { + char filename[100]; #ifdef PARALLEL - sprintf (filename, "netgen.prof.%d", id); + sprintf (filename, "netgen.prof.%d", id); #else - sprintf (filename, "netgen.prof"); + sprintf (filename, "netgen.prof"); #endif - - FILE *prof = fopen(filename,"w"); - Print (prof); - fclose(prof); - */ + + printf ("write profile to file %s\n", filename); + FILE *prof = fopen(filename,"w"); + Print (prof); + fclose(prof); + } } diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 4292b788..4ec3dab1 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -2643,6 +2643,10 @@ namespace netgen PrintMessage(1,"Mesh bisection"); PushStatus("Mesh bisection"); + static int timer = NgProfiler::CreateTimer ("Bisect"); + NgProfiler::RegionTimer reg1 (timer); + + static int localizetimer = NgProfiler::CreateTimer("localize edgepoints"); NgProfiler::RegionTimer * loct = new NgProfiler::RegionTimer(localizetimer); @@ -3524,6 +3528,7 @@ namespace netgen } while (hangingvol || hangingsurf || hangingedge); + /* if (printmessage_importance>0) { ostringstream strstr; @@ -3536,8 +3541,16 @@ namespace netgen } strstr << mesh.GetNP() << " points"; PrintMessage(4,strstr.str()); - } + */ + PrintMessage (4, mtets.Size(), " tets"); + PrintMessage (4, mtris.Size(), " trigs"); + if (mprisms.Size()) + { + PrintMessage (4, mprisms.Size(), " prisms"); + PrintMessage (4, mquads.Size(), " quads"); + } + PrintMessage (4, mesh.GetNP(), " points"); } @@ -3859,9 +3872,6 @@ namespace netgen // Check/Repair - //cout << "Hallo Welt" << endl; - //getchar(); - static bool repaired_once; if(mesh.mglevels == 1) repaired_once = false; diff --git a/libsrc/meshing/classifyhpel.hpp b/libsrc/meshing/classifyhpel.hpp index 69656ffc..c1555371 100644 --- a/libsrc/meshing/classifyhpel.hpp +++ b/libsrc/meshing/classifyhpel.hpp @@ -455,7 +455,7 @@ HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE & edg point_sing[p[j]-1] = 1; } - const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (PRISM); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (PRISM); for(int k=0;k<9;k++) { INDEX_2 i2 = INDEX_2 :: Sort(el.PNum(p[eledges[k][0]-1]),el.PNum(p[eledges[k][1]-1])); @@ -463,7 +463,7 @@ HPREF_ELEMENT_TYPE ClassifyPrism(HPRefElement & el, INDEX_2_HASHTABLE & edg else edge_sing[k] = face_edges.Used(i2); } - const ELEMENT_FACE * elfaces = MeshTopology::GetFaces (PRISM); + const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (PRISM); for (int k=0;k<5;k++) { INDEX_3 i3; @@ -710,7 +710,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge *testout << endl; */ } - const ELEMENT_EDGE * eledges = MeshTopology::GetEdges(TRIG); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1(TRIG); if(dim==3) { @@ -1499,8 +1499,8 @@ HPREF_ELEMENT_TYPE ClassifyHex(HPRefElement & el, INDEX_2_HASHTABLE & edges // indices of bot,top-faces combinations int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; int p[8]; - const ELEMENT_FACE * elfaces = MeshTopology::GetFaces (HEX); - const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (HEX); + const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (HEX); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (HEX); for(int m=0;m<6 && type == HP_NONE;m++) for(int j=0;j<4 && type == HP_NONE;j++) @@ -1643,8 +1643,8 @@ HPREF_ELEMENT_TYPE ClassifyPyramid(HPRefElement & el, INDEX_2_HASHTABLE & e // indices of bot,top-faces combinations // int index[6][2] = {{0,1},{1,0},{2,4},{4,2},{3,5},{5,3}}; - const ELEMENT_FACE * elfaces = MeshTopology::GetFaces (PYRAMID); - const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (PYRAMID); + const ELEMENT_FACE * elfaces = MeshTopology::GetFaces1 (PYRAMID); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (PYRAMID); int point_sing[5]={0,0,0,0,0}; int face_sing[5] = {0,0,0,0,0}; diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp index 412ec8b4..3282613c 100644 --- a/libsrc/meshing/clusters.cpp +++ b/libsrc/meshing/clusters.cpp @@ -5,263 +5,263 @@ namespace netgen { -AnisotropicClusters :: AnisotropicClusters (const Mesh & amesh) - : mesh(amesh) -{ - ; -} + AnisotropicClusters :: AnisotropicClusters (const Mesh & amesh) + : mesh(amesh) + { + ; + } -AnisotropicClusters :: ~AnisotropicClusters () -{ - ; -} + AnisotropicClusters :: ~AnisotropicClusters () + { + ; + } -void AnisotropicClusters :: Update() -{ - int i, j, k; + void AnisotropicClusters :: Update() + { + int i, j, k; - const MeshTopology & top = mesh.GetTopology(); + const MeshTopology & top = mesh.GetTopology(); - bool hasedges = top.HasEdges(); - bool hasfaces = top.HasFaces(); + bool hasedges = top.HasEdges(); + bool hasfaces = top.HasFaces(); - if (!hasedges || !hasfaces) return; + if (!hasedges || !hasfaces) return; - PrintMessage (3, "Update Clusters"); + PrintMessage (3, "Update Clusters"); - nv = mesh.GetNV(); - ned = top.GetNEdges(); - nfa = top.GetNFaces(); - ne = mesh.GetNE(); - int nse = mesh.GetNSE(); + nv = mesh.GetNV(); + ned = top.GetNEdges(); + nfa = top.GetNFaces(); + ne = mesh.GetNE(); + int nse = mesh.GetNSE(); - cluster_reps.SetSize (nv+ned+nfa+ne); + cluster_reps.SetSize (nv+ned+nfa+ne); - Array nnums, ednums, fanums; - int changed; + Array nnums, ednums, fanums; + int changed; - for (i = 1; i <= cluster_reps.Size(); i++) - cluster_reps.Elem(i) = -1; + for (i = 1; i <= cluster_reps.Size(); i++) + cluster_reps.Elem(i) = -1; - for (i = 1; i <= ne; i++) - { - const Element & el = mesh.VolumeElement(i); - ELEMENT_TYPE typ = el.GetType(); + for (i = 1; i <= ne; i++) + { + const Element & el = mesh.VolumeElement(i); + ELEMENT_TYPE typ = el.GetType(); - top.GetElementEdges (i, ednums); - top.GetElementFaces (i, fanums); + top.GetElementEdges (i, ednums); + top.GetElementFaces (i, fanums); - int elnv = top.GetNVertices (typ); - int elned = ednums.Size(); - int elnfa = fanums.Size(); + int elnv = top.GetNVertices (typ); + int elned = ednums.Size(); + int elnfa = fanums.Size(); - nnums.SetSize(elnv+elned+elnfa+1); - for (j = 1; j <= elnv; j++) - nnums.Elem(j) = el.PNum(j); - for (j = 1; j <= elned; j++) - nnums.Elem(elnv+j) = nv+ednums.Elem(j); - for (j = 1; j <= elnfa; j++) - nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j); - nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i; + nnums.SetSize(elnv+elned+elnfa+1); + for (j = 1; j <= elnv; j++) + nnums.Elem(j) = el.PNum(j); + for (j = 1; j <= elned; j++) + nnums.Elem(elnv+j) = nv+ednums.Elem(j); + for (j = 1; j <= elnfa; j++) + nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j); + nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i; - for (j = 0; j < nnums.Size(); j++) - cluster_reps.Elem(nnums[j]) = nnums[j]; - } + for (j = 0; j < nnums.Size(); j++) + cluster_reps.Elem(nnums[j]) = nnums[j]; + } - for (i = 1; i <= nse; i++) - { - const Element2d & el = mesh.SurfaceElement(i); - ELEMENT_TYPE typ = el.GetType(); + for (i = 1; i <= nse; i++) + { + const Element2d & el = mesh.SurfaceElement(i); + ELEMENT_TYPE typ = el.GetType(); - top.GetSurfaceElementEdges (i, ednums); - int fanum = top.GetSurfaceElementFace (i); + top.GetSurfaceElementEdges (i, ednums); + int fanum = top.GetSurfaceElementFace (i); - int elnv = top.GetNVertices (typ); - int elned = ednums.Size(); + int elnv = top.GetNVertices (typ); + int elned = ednums.Size(); - nnums.SetSize(elnv+elned+1); - for (j = 1; j <= elnv; j++) - nnums.Elem(j) = el.PNum(j); - for (j = 1; j <= elned; j++) - nnums.Elem(elnv+j) = nv+ednums.Elem(j); - nnums.Elem(elnv+elned+1) = fanum; + nnums.SetSize(elnv+elned+1); + for (j = 1; j <= elnv; j++) + nnums.Elem(j) = el.PNum(j); + for (j = 1; j <= elned; j++) + nnums.Elem(elnv+j) = nv+ednums.Elem(j); + nnums.Elem(elnv+elned+1) = fanum; - for (j = 0; j < nnums.Size(); j++) - cluster_reps.Elem(nnums[j]) = nnums[j]; - } + for (j = 0; j < nnums.Size(); j++) + cluster_reps.Elem(nnums[j]) = nnums[j]; + } - static const int hex_cluster[] = - { - 1, 2, 3, 4, 1, 2, 3, 4, - 5, 6, 7, 8, 5, 6, 7, 8, 1, 2, 3, 4, - 9, 9, 5, 8, 6, 7, - 9 - }; + static const int hex_cluster[] = + { + 1, 2, 3, 4, 1, 2, 3, 4, + 5, 6, 7, 8, 5, 6, 7, 8, 1, 2, 3, 4, + 9, 9, 5, 8, 6, 7, + 9 + }; - static const int prism_cluster[] = - { - 1, 2, 3, 1, 2, 3, - 4, 5, 6, 4, 5, 6, 3, 1, 2, - 7, 7, 4, 5, 6, - 7 - }; + static const int prism_cluster[] = + { + 1, 2, 3, 1, 2, 3, + 4, 5, 6, 4, 5, 6, 3, 1, 2, + 7, 7, 4, 5, 6, + 7 + }; - static const int pyramid_cluster[] = - { - 1, 2, 2, 1, 3, - 4, 2, 1, 4, 6, 5, 5, 6, - 7, 5, 7, 6, 4, - 7 - }; - static const int tet_cluster14[] = - { 1, 2, 3, 1, 1, 4, 5, 4, 5, 6, 7, 5, 4, 7, 7 }; + static const int pyramid_cluster[] = + { + 1, 2, 2, 1, 3, + 4, 2, 1, 4, 6, 5, 5, 6, + 7, 5, 7, 6, 4, + 7 + }; + static const int tet_cluster14[] = + { 1, 2, 3, 1, 1, 4, 5, 4, 5, 6, 7, 5, 4, 7, 7 }; - static const int tet_cluster12[] = - { 1, 1, 2, 3, 4, 4, 5, 1, 6, 6, 7, 7, 4, 6, 7 }; + static const int tet_cluster12[] = + { 1, 1, 2, 3, 4, 4, 5, 1, 6, 6, 7, 7, 4, 6, 7 }; - static const int tet_cluster13[] = - { 1, 2, 1, 3, 4, 6, 4, 5, 1, 5, 7, 4, 7, 5, 7 }; + static const int tet_cluster13[] = + { 1, 2, 1, 3, 4, 6, 4, 5, 1, 5, 7, 4, 7, 5, 7 }; - static const int tet_cluster23[] = - { 2, 1, 1, 3, 6, 5, 5, 4, 4, 1, 5, 7, 7, 4, 7 }; + static const int tet_cluster23[] = + { 2, 1, 1, 3, 6, 5, 5, 4, 4, 1, 5, 7, 7, 4, 7 }; - static const int tet_cluster24[] = - { 2, 1, 3, 1, 4, 1, 5, 4, 6, 5, 5, 7, 4, 7, 7 }; + static const int tet_cluster24[] = + { 2, 1, 3, 1, 4, 1, 5, 4, 6, 5, 5, 7, 4, 7, 7 }; - static const int tet_cluster34[] = - { 2, 3, 1, 1, 4, 5, 1, 6, 4, 5, 5, 4, 7, 7, 7 }; + static const int tet_cluster34[] = + { 2, 3, 1, 1, 4, 5, 1, 6, 4, 5, 5, 4, 7, 7, 7 }; - int cnt = 0; + int cnt = 0; - do - { + do + { - cnt++; - changed = 0; + cnt++; + changed = 0; - for (i = 1; i <= ne; i++) - { - const Element & el = mesh.VolumeElement(i); - ELEMENT_TYPE typ = el.GetType(); + for (i = 1; i <= ne; i++) + { + const Element & el = mesh.VolumeElement(i); + ELEMENT_TYPE typ = el.GetType(); - top.GetElementEdges (i, ednums); - top.GetElementFaces (i, fanums); + top.GetElementEdges (i, ednums); + top.GetElementFaces (i, fanums); - int elnv = top.GetNVertices (typ); - int elned = ednums.Size(); - int elnfa = fanums.Size(); + int elnv = top.GetNVertices (typ); + int elned = ednums.Size(); + int elnfa = fanums.Size(); - nnums.SetSize(elnv+elned+elnfa+1); - for (j = 1; j <= elnv; j++) - nnums.Elem(j) = el.PNum(j); - for (j = 1; j <= elned; j++) - nnums.Elem(elnv+j) = nv+ednums.Elem(j); - for (j = 1; j <= elnfa; j++) - nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j); - nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i; + nnums.SetSize(elnv+elned+elnfa+1); + for (j = 1; j <= elnv; j++) + nnums.Elem(j) = el.PNum(j); + for (j = 1; j <= elned; j++) + nnums.Elem(elnv+j) = nv+ednums.Elem(j); + for (j = 1; j <= elnfa; j++) + nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j); + nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i; - const int * clustertab = NULL; - switch (typ) - { - case PRISM: - case PRISM12: - clustertab = prism_cluster; - break; - case HEX: - clustertab = hex_cluster; - break; - case PYRAMID: - clustertab = pyramid_cluster; - break; - case TET: - case TET10: - if (cluster_reps.Get(el.PNum(1)) == - cluster_reps.Get(el.PNum(2))) - clustertab = tet_cluster12; - else if (cluster_reps.Get(el.PNum(1)) == - cluster_reps.Get(el.PNum(3))) - clustertab = tet_cluster13; - else if (cluster_reps.Get(el.PNum(1)) == - cluster_reps.Get(el.PNum(4))) - clustertab = tet_cluster14; - else if (cluster_reps.Get(el.PNum(2)) == - cluster_reps.Get(el.PNum(3))) - clustertab = tet_cluster23; - else if (cluster_reps.Get(el.PNum(2)) == - cluster_reps.Get(el.PNum(4))) - clustertab = tet_cluster24; - else if (cluster_reps.Get(el.PNum(3)) == - cluster_reps.Get(el.PNum(4))) - clustertab = tet_cluster34; + const int * clustertab = NULL; + switch (typ) + { + case PRISM: + case PRISM12: + clustertab = prism_cluster; + break; + case HEX: + clustertab = hex_cluster; + break; + case PYRAMID: + clustertab = pyramid_cluster; + break; + case TET: + case TET10: + if (cluster_reps.Get(el.PNum(1)) == + cluster_reps.Get(el.PNum(2))) + clustertab = tet_cluster12; + else if (cluster_reps.Get(el.PNum(1)) == + cluster_reps.Get(el.PNum(3))) + clustertab = tet_cluster13; + else if (cluster_reps.Get(el.PNum(1)) == + cluster_reps.Get(el.PNum(4))) + clustertab = tet_cluster14; + else if (cluster_reps.Get(el.PNum(2)) == + cluster_reps.Get(el.PNum(3))) + clustertab = tet_cluster23; + else if (cluster_reps.Get(el.PNum(2)) == + cluster_reps.Get(el.PNum(4))) + clustertab = tet_cluster24; + else if (cluster_reps.Get(el.PNum(3)) == + cluster_reps.Get(el.PNum(4))) + clustertab = tet_cluster34; - else + else + clustertab = NULL; + break; + default: clustertab = NULL; - break; - default: - clustertab = NULL; - } + } - if (clustertab) - for (j = 0; j < nnums.Size(); j++) - for (k = 0; k < j; k++) - if (clustertab[j] == clustertab[k]) - { - int jj = nnums[j]; - int kk = nnums[k]; - if (cluster_reps.Get(jj) < cluster_reps.Get(kk)) - { - cluster_reps.Elem(kk) = cluster_reps.Get(jj); - changed = 1; - } - else if (cluster_reps.Get(kk) < cluster_reps.Get(jj)) - { - cluster_reps.Elem(jj) = cluster_reps.Get(kk); - changed = 1; - } - } + if (clustertab) + for (j = 0; j < nnums.Size(); j++) + for (k = 0; k < j; k++) + if (clustertab[j] == clustertab[k]) + { + int jj = nnums[j]; + int kk = nnums[k]; + if (cluster_reps.Get(jj) < cluster_reps.Get(kk)) + { + cluster_reps.Elem(kk) = cluster_reps.Get(jj); + changed = 1; + } + else if (cluster_reps.Get(kk) < cluster_reps.Get(jj)) + { + cluster_reps.Elem(jj) = cluster_reps.Get(kk); + changed = 1; + } + } - /* - if (clustertab) - { + /* + if (clustertab) + { if (typ == PYRAMID) - (*testout) << "pyramid"; + (*testout) << "pyramid"; else if (typ == PRISM || typ == PRISM12) - (*testout) << "prism"; + (*testout) << "prism"; else if (typ == TET || typ == TET10) - (*testout) << "tet"; + (*testout) << "tet"; else - (*testout) << "unknown type" << endl; + (*testout) << "unknown type" << endl; (*testout) << ", nnums = "; for (j = 0; j < nnums.Size(); j++) - (*testout) << "node " << j << " = " << nnums[j] << ", rep = " - << cluster_reps.Get(nnums[j]) << endl; - } - */ - } - } - while (changed); + (*testout) << "node " << j << " = " << nnums[j] << ", rep = " + << cluster_reps.Get(nnums[j]) << endl; + } + */ + } + } + while (changed); - /* - (*testout) << "cluster reps:" << endl; - for (i = 1; i <= cluster_reps.Size(); i++) - { - (*testout) << i << ": "; - if (i <= nv) - (*testout) << "v" << i << " "; - else if (i <= nv+ned) - (*testout) << "e" << i-nv << " "; - else if (i <= nv+ned+nfa) - (*testout) << "f" << i-nv-ned << " "; - else - (*testout) << "c" << i-nv-ned-nfa << " "; - (*testout) << cluster_reps.Get(i) << endl; - } - */ -} + /* + (*testout) << "cluster reps:" << endl; + for (i = 1; i <= cluster_reps.Size(); i++) + { + (*testout) << i << ": "; + if (i <= nv) + (*testout) << "v" << i << " "; + else if (i <= nv+ned) + (*testout) << "e" << i-nv << " "; + else if (i <= nv+ned+nfa) + (*testout) << "f" << i-nv-ned << " "; + else + (*testout) << "c" << i-nv-ned-nfa << " "; + (*testout) << cluster_reps.Get(i) << endl; + } + */ + } } diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index cdc05901..51c29fb9 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -545,7 +545,7 @@ namespace netgen top.GetSurfaceElementEdges (i+1, edgenrs); for (int j = 0; j < edgenrs.Size(); j++) edgenrs[j]--; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (el.GetType()); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (el.GetType()); for (int i2 = 0; i2 < edgenrs.Size(); i2++) { @@ -1371,7 +1371,7 @@ namespace netgen for (int j = 0; j < 3; j++) shapes(j) = lami[j] * lami[j]; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG); for (int j = 0; j < 3; j++) { double wi = edgeweight[info.edgenrs[j]]; @@ -1394,7 +1394,7 @@ namespace netgen if (info.order == 1) return; int ii = 3; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG); for (int i = 0; i < 3; i++) { @@ -1465,7 +1465,7 @@ namespace netgen }; int ii = 4; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (QUAD); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (QUAD); for (int i = 0; i < 4; i++) { @@ -1525,7 +1525,7 @@ namespace netgen dshapes(j,1) = 2 * lami[j] * dlami[j][1]; } - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG); for (int j = 0; j < 3; j++) { double wi = edgeweight[info.edgenrs[j]]; @@ -1573,7 +1573,7 @@ namespace netgen lami[2] = 1-xi(0)-xi(1); int ii = 3; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TRIG); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TRIG); for (int i = 0; i < 3; i++) { @@ -1707,7 +1707,7 @@ namespace netgen ArrayMem hshapes(order+1), hdshapes(order+1); int ii = 4; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (QUAD); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (QUAD); for (int i = 0; i < 4; i++) { @@ -1987,7 +1987,7 @@ namespace netgen for (int j = 0; j < 4; j++) shapes(j) = lami[j] * lami[j]; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET); for (int j = 0; j < 6; j++) { double wi = edgeweight[info.edgenrs[j]]; @@ -2013,7 +2013,7 @@ namespace netgen if (info.order == 1) return; int ii = 4; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET); for (int i = 0; i < 6; i++) { int eorder = edgeorder[info.edgenrs[i]]; @@ -2026,7 +2026,7 @@ namespace netgen ii += eorder-1; } } - const ELEMENT_FACE * faces = MeshTopology::GetFaces (TET); + const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (TET); for (int i = 0; i < 4; i++) { int forder = faceorder[info.facenrs[i]]; @@ -2091,7 +2091,7 @@ namespace netgen int ii = 6; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (PRISM); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (PRISM); for (int i = 0; i < 6; i++) // horizontal edges { int eorder = edgeorder[info.edgenrs[i]]; @@ -2131,7 +2131,7 @@ namespace netgen } // FACE SHAPES - const ELEMENT_FACE * faces = MeshTopology::GetFaces (PRISM); + const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (PRISM); for (int i = 0; i < 2; i++) { int forder = faceorder[info.facenrs[i]]; @@ -2170,7 +2170,7 @@ namespace netgen if (info.order == 1) return; int ii = 5; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (PYRAMID); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (PYRAMID); for (int i = 0; i < 4; i++) // horizontal edges { int eorder = edgeorder[info.edgenrs[i]]; @@ -2245,7 +2245,7 @@ namespace netgen dshapes(j,2) = 2 * lami[j] * dlami[j][2]; } - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET); for (int j = 0; j < 6; j++) { double wi = edgeweight[info.edgenrs[j]]; @@ -2283,7 +2283,7 @@ namespace netgen double lami[] = { xi(0), xi(1), xi(2), 1-xi(0)-xi(1)-xi(2) }; int ii = 4; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (TET); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (TET); for (int i = 0; i < 6; i++) { int eorder = edgeorder[info.edgenrs[i]]; @@ -2314,7 +2314,7 @@ namespace netgen } } - const ELEMENT_FACE * faces = MeshTopology::GetFaces (TET); + const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (TET); for (int i = 0; i < 4; i++) { int forder = faceorder[info.facenrs[i]]; @@ -2428,7 +2428,7 @@ namespace netgen if (info.order == 1) return; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (PRISM); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (PRISM); for (int i = 0; i < 6; i++) // horizontal edges { int order = edgeorder[info.edgenrs[i]]; @@ -2502,7 +2502,7 @@ namespace netgen if (info.order == 2) return; // FACE SHAPES - const ELEMENT_FACE * faces = MeshTopology::GetFaces (PRISM); + const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (PRISM); for (int i = 0; i < 2; i++) { int forder = faceorder[info.facenrs[i]]; diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index ca90f129..9c8155ad 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -991,7 +991,7 @@ namespace netgen int npi = mesh.AddPoint (center); - const ELEMENT_FACE * faces = MeshTopology::GetFaces (HEX); + const ELEMENT_FACE * faces = MeshTopology::GetFaces1 (HEX); for (int j = 0; j < 6; j++) { @@ -1448,7 +1448,7 @@ namespace netgen { Element el = mesh[i] ; HPRefElement & hpel = hpelements[mesh[i].hp_elnr]; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (mesh[i].GetType()); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType()); double dist[3] = {0,0,0}; int ord_dir[3] = {0,0,0}; int edge_dir[12] = {0,0,0,0,0,0,0,0,0,0,0,0}; @@ -1520,7 +1520,7 @@ namespace netgen { Element2d el = mesh[i] ; HPRefElement & hpel = hpelements[mesh[i].hp_elnr]; - const ELEMENT_EDGE * edges = MeshTopology::GetEdges (mesh[i].GetType()); + const ELEMENT_EDGE * edges = MeshTopology::GetEdges1 (mesh[i].GetType()); double dist[3] = {0,0,0}; int ord_dir[3] = {0,0,0}; int edge_dir[4] = {0,0,0,0} ; @@ -1637,7 +1637,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS for (int i = 1; i <= mesh.GetNE(); i++) { const Element & el = mesh.VolumeElement(i); - const ELEMENT_EDGE * eledges = MeshTopology::GetEdges (el.GetType()); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges1 (el.GetType()); int nedges = MeshTopology::GetNEdges (el.GetType()); for (int j = 0; j < nedges; j++) for (int k = 0; k < nedges; k++) diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 65750171..d5b0dcf0 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -2,812 +2,1155 @@ #include "meshing.hpp" +#define DEBUG namespace netgen { -MeshTopology :: MeshTopology (const Mesh & amesh) - : mesh(amesh) -{ - buildedges = 1; - buildfaces = 1; - vert2element = 0; - vert2surfelement = 0; - vert2segment = 0; - timestamp = -1; + MeshTopology :: MeshTopology (const Mesh & amesh) + : mesh(amesh) + { + buildedges = 1; + buildfaces = 1; + vert2element = 0; + vert2surfelement = 0; + vert2segment = 0; + timestamp = -1; - edge2vert.SetName ("edge2vert"); - face2vert.SetName ("face2vert"); - edges.SetName ("el2edge"); - faces.SetName ("el2face"); - surfedges.SetName ("surfel2edge"); - segedges.SetName ("segment2edge"); - surffaces.SetName ("surfel2face"); - surf2volelement.SetName ("surfel2el"); - face2surfel.SetName ("face2surfel"); -} + edge2vert.SetName ("edge2vert"); + face2vert.SetName ("face2vert"); + edges.SetName ("el2edge"); + faces.SetName ("el2face"); + surfedges.SetName ("surfel2edge"); + segedges.SetName ("segment2edge"); + surffaces.SetName ("surfel2face"); + surf2volelement.SetName ("surfel2el"); + face2surfel.SetName ("face2surfel"); + } -MeshTopology :: ~MeshTopology () -{ - delete vert2element; - delete vert2surfelement; - delete vert2segment; -} + MeshTopology :: ~MeshTopology () + { + delete vert2element; + delete vert2surfelement; + delete vert2segment; + } -void MeshTopology :: Update() -{ - static int timer = NgProfiler::CreateTimer ("topology"); - NgProfiler::RegionTimer reg (timer); + void MeshTopology :: Update() + { + static int timer = NgProfiler::CreateTimer ("topology"); + NgProfiler::RegionTimer reg (timer); #ifdef PARALLEL - ParallelMeshTopology & paralleltop = mesh.GetParallelTopology(); + ParallelMeshTopology & paralleltop = mesh.GetParallelTopology(); - bool isparallel = 0; + bool isparallel = 0; #endif - if (timestamp > mesh.GetTimeStamp()) return; + if (timestamp > mesh.GetTimeStamp()) return; - int ne = mesh.GetNE(); - int nse = mesh.GetNSE(); - int nseg = mesh.GetNSeg(); - int np = mesh.GetNP(); - int nv = mesh.GetNV(); - int nfa = 0; - int ned = edge2vert.Size(); + int ne = mesh.GetNE(); + int nse = mesh.GetNSE(); + int nseg = mesh.GetNSeg(); + int np = mesh.GetNP(); + int nv = mesh.GetNV(); + int nfa = 0; + int ned = edge2vert.Size(); - PrintMessage (3, "Update mesh topology"); + PrintMessage (3, "Update mesh topology"); - (*testout) << " UPDATE MESH TOPOLOGY " << endl; - (*testout) << "ne = " << ne << endl; - (*testout) << "nse = " << nse << endl; - (*testout) << "nseg = " << nseg << endl; - (*testout) << "np = " << np << endl; - (*testout) << "nv = " << nv << endl; + (*testout) << " UPDATE MESH TOPOLOGY " << endl; + (*testout) << "ne = " << ne << endl; + (*testout) << "nse = " << nse << endl; + (*testout) << "nseg = " << nseg << endl; + (*testout) << "np = " << np << endl; + (*testout) << "nv = " << nv << endl; - delete vert2element; - delete vert2surfelement; - delete vert2segment; + delete vert2element; + delete vert2surfelement; + delete vert2segment; - Array cnt(nv); - Array vnums; + Array cnt(nv); + Array vnums; - /* - generate: - vertex to element - vertex to surface element - vertex to segment - */ + /* + generate: + vertex to element + vertex to surface element + vertex to segment + */ - cnt = 0; - for (ElementIndex ei = 0; ei < ne; ei++) - { - const Element & el = mesh[ei]; - int nelv = el.GetNV(); - for (int j = 0; j < nelv; j++) - cnt[el[j]]++; - } + cnt = 0; + for (ElementIndex ei = 0; ei < ne; ei++) + { + const Element & el = mesh[ei]; + int nelv = el.GetNV(); + for (int j = 0; j < nelv; j++) + cnt[el[j]]++; + } - vert2element = new TABLE (cnt); - for (ElementIndex ei = 0; ei < ne; ei++) - { - const Element & el = mesh[ei]; - int nelv = el.GetNV(); - for (int j = 0; j < nelv; j++) - vert2element->AddSave (el[j], ei+1); - } + vert2element = new TABLE (cnt); + for (ElementIndex ei = 0; ei < ne; ei++) + { + const Element & el = mesh[ei]; + int nelv = el.GetNV(); + for (int j = 0; j < nelv; j++) + vert2element->AddSave (el[j], ei+1); + } - cnt = 0; - for (SurfaceElementIndex sei = 0; sei < nse; sei++) - { - const Element2d & el = mesh[sei]; - int nelv = el.GetNV(); - for (int j = 0; j < nelv; j++) - cnt[el[j]]++; - } + cnt = 0; + for (SurfaceElementIndex sei = 0; sei < nse; sei++) + { + const Element2d & el = mesh[sei]; + int nelv = el.GetNV(); + for (int j = 0; j < nelv; j++) + cnt[el[j]]++; + } - vert2surfelement = new TABLE (cnt); - for (SurfaceElementIndex sei = 0; sei < nse; sei++) - { - const Element2d & el = mesh[sei]; - int nelv = el.GetNV(); - for (int j = 0; j < nelv; j++) - vert2surfelement->AddSave (el[j], sei+1); - } + vert2surfelement = new TABLE (cnt); + for (SurfaceElementIndex sei = 0; sei < nse; sei++) + { + const Element2d & el = mesh[sei]; + int nelv = el.GetNV(); + for (int j = 0; j < nelv; j++) + vert2surfelement->AddSave (el[j], sei+1); + } - cnt = 0; - for (int i = 1; i <= nseg; i++) - { - const Segment & seg = mesh.LineSegment(i); - cnt[seg[0]]++; - cnt[seg[1]]++; - } + cnt = 0; + for (int i = 1; i <= nseg; i++) + { + const Segment & seg = mesh.LineSegment(i); + cnt[seg[0]]++; + cnt[seg[1]]++; + } - vert2segment = new TABLE (cnt); - for (int i = 1; i <= nseg; i++) - { - const Segment & seg = mesh.LineSegment(i); - vert2segment->AddSave (seg[0], i); - vert2segment->AddSave (seg[1], i); - } + vert2segment = new TABLE (cnt); + for (int i = 1; i <= nseg; i++) + { + const Segment & seg = mesh.LineSegment(i); + vert2segment->AddSave (seg[0], i); + vert2segment->AddSave (seg[1], i); + } - if (buildedges) - { - static int timer1 = NgProfiler::CreateTimer ("topology::buildedges"); - NgProfiler::RegionTimer reg1 (timer1); + - PrintMessage (5, "Update edges "); + if (buildedges) + { + static int timer1 = NgProfiler::CreateTimer ("topology::buildedges"); + NgProfiler::RegionTimer reg1 (timer1); + + PrintMessage (5, "Update edges "); - edges.SetSize(ne); - surfedges.SetSize(nse); - segedges.SetSize(nseg); + edges.SetSize(ne); + surfedges.SetSize(nse); + segedges.SetSize(nseg); - for (int i = 0; i < ne; i++) - for (int j = 0; j < 12; j++) - edges[i][j] = 0; - for (int i = 0; i < nse; i++) - for (int j = 0; j < 4; j++) - surfedges[i][j] = 0; + for (int i = 0; i < ne; i++) + for (int j = 0; j < 12; j++) + edges[i][j] = 0; + for (int i = 0; i < nse; i++) + for (int j = 0; j < 4; j++) + surfedges[i][j] = 0; - // keep existing edges - cnt = 0; - for (int i = 0; i < edge2vert.Size(); i++) - cnt[edge2vert[i][0]]++; - TABLE vert2edge (cnt); - for (int i = 0; i < edge2vert.Size(); i++) - vert2edge.AddSave (edge2vert[i][0], i+1); + // keep existing edges + cnt = 0; + for (int i = 0; i < edge2vert.Size(); i++) + cnt[edge2vert[i][0]]++; + TABLE vert2edge (cnt); + for (int i = 0; i < edge2vert.Size(); i++) + vert2edge.AddSave (edge2vert[i][0], i+1); - // ensure all coarse grid and intermediate level edges - cnt = 0; - for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) - { - int pa[2]; - pa[0] = mesh.mlbetweennodes[i].I1(); - pa[1] = mesh.mlbetweennodes[i].I2(); - if (pa[0] > pa[1]) Swap (pa[0], pa[1]); - if (pa[0] > 0) - cnt.Elem(pa[0])++; - } - TABLE vert2vertcoarse (cnt); - for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) - { - int pa[2]; - pa[0] = mesh.mlbetweennodes[i].I1(); - pa[1] = mesh.mlbetweennodes[i].I2(); - if (pa[0] > pa[1]) swap (pa[0], pa[1]); - if (pa[0] > 0) - vert2vertcoarse.AddSave1 (pa[0], pa[1]); - } - - - Array edgenr(nv); - Array edgeflag(nv); - edgeflag = 0; - - ned = edge2vert.Size(); - Array missing; - - for (int i = 1; i <= nv; i++) - { - for (int j = 1; j <= vert2edge.EntrySize(i); j++) - { - int ednr = vert2edge.Get(i,j); - int i2 = edge2vert.Get(ednr)[1]; - edgeflag[i2] = i; - edgenr[i2] = ednr; - } - for (int j = 1; j <= vert2vertcoarse.EntrySize(i); j++) - { - int v2 = vert2vertcoarse.Get(i,j); - if (edgeflag[v2] < i) - { - ned++; - edgenr[v2] = ned; - edgeflag[v2] = i; - missing.Append (INDEX_3(i,v2,ned)); - } - } - - for (int j = 1; j <= vert2element->EntrySize(i); j++) - { - int elnr = vert2element->Get(i,j); - const Element & el = mesh.VolumeElement (elnr); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el.PNum(eledges[k][0]), - el.PNum(eledges[k][1])); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) - continue; - - if (edgeflag[edge.I2()] < i) - { - ned++; - edgenr[edge.I2()] = ned; - edgeflag[edge.I2()] = i; - } - - int edgenum = edgenr[edge.I2()]; - if (edgedir) edgenum *= -1; - edges.Elem(elnr)[k] = edgenum; - } - } - - for (int j = 1; j <= vert2surfelement->EntrySize(i); j++) - { - int elnr = vert2surfelement->Get(i,j); - const Element2d & el = mesh.SurfaceElement (elnr); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el.PNum(eledges[k][0]), - el.PNum(eledges[k][1])); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) - continue; - - if (edgeflag[edge.I2()] < i) - { - ned++; - edgenr[edge.I2()] = ned; - edgeflag[edge.I2()] = i; - } - - int edgenum = edgenr[edge.I2()]; - if (edgedir) edgenum *= -1; - surfedges.Elem(elnr)[k] = edgenum; - } - } - - for (int j = 1; j <= vert2segment->EntrySize(i); j++) - { - int elnr = vert2segment->Get(i,j); - const Segment & el = mesh.LineSegment (elnr); - - INDEX_2 edge(el[0], el[1]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) - continue; - - if (edgeflag[edge.I2()] < i) - { - ned++; - edgenr[edge.I2()] = ned; - edgeflag[edge.I2()] = i; - } - int edgenum = edgenr[edge.I2()]; - - if (edgedir) edgenum *= -1; - segedges.Elem(elnr) = edgenum; - } - } - - - edge2vert.SetSize (ned); - for (int i = 1; i <= ne; i++) - { - const Element & el = mesh.VolumeElement (i); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el.PNum(eledges[k][0]), - el.PNum(eledges[k][1])); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - int edgenum = abs (edges.Elem(i)[k]); - - edge2vert.Elem(edgenum)[0] = edge.I1(); - edge2vert.Elem(edgenum)[1] = edge.I2(); - } - } - for (int i = 1; i <= nse; i++) - { - const Element2d & el = mesh.SurfaceElement (i); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el.PNum(eledges[k][0]), - el.PNum(eledges[k][1])); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - int edgenum = abs (surfedges.Elem(i)[k]); - - edge2vert.Elem(edgenum)[0] = edge.I1(); - edge2vert.Elem(edgenum)[1] = edge.I2(); - } - } - - for (int i = 1; i <= nseg; i++) - { - const Segment & el = mesh.LineSegment (i); - - INDEX_2 edge(el[0], el[1]); - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - int edgenum = abs (segedges.Elem(i)); - - edge2vert.Elem(edgenum)[0] = edge.I1(); - edge2vert.Elem(edgenum)[1] = edge.I2(); - } - - for (int i = 1; i <= missing.Size(); i++) - { - INDEX_3 i3 = missing.Get(i); - edge2vert.Elem(i3.I3())[0] = i3.I1(); - edge2vert.Elem(i3.I3())[1] = i3.I2(); - } - - - /* - (*testout) << "edge table:" << endl; - (*testout) << "edge2vert:" << endl; - for (int i = 1; i <= edge2vert.Size(); i++) - (*testout) << "edge " << i << ", v1,2 = " << edge2vert.Elem(i)[0] << ", " << edge2vert.Elem(i)[1] << endl; - (*testout) << "surfedges:" << endl; - for (int i = 1; i <= surfedges.Size(); i++) - (*testout) << "el " << i << ", edges = " - << surfedges.Elem(i)[0] << ", " - << surfedges.Elem(i)[1] << ", " - << surfedges.Elem(i)[2] << endl; - */ - - - } - - - // cout << "build edges done" << endl; - - // generate faces - if (buildfaces) // && mesh.GetDimension() == 3) - { - int i, j; - - static int timer2 = NgProfiler::CreateTimer ("topology::buildfaces"); - NgProfiler::RegionTimer reg2 (timer2); - - PrintMessage (5, "Update faces "); - - faces.SetSize(ne); - surffaces.SetSize(nse); - - // face2vert.SetSize(0); // keep old faces - nfa = face2vert.Size(); - // INDEX_3_HASHTABLE vert2face(ne+nse+1); - INDEX_3_CLOSED_HASHTABLE vert2face(8*ne+2*nse+nfa+2); - - for (i = 1; i <= face2vert.Size(); i++) - { - INDEX_3 f; - f.I1() = face2vert.Get(i)[0]; - f.I2() = face2vert.Get(i)[1]; - f.I3() = face2vert.Get(i)[2]; - vert2face.Set (f, i); - } - - for (i = 1; i <= ne; i++) - { - const Element & el = mesh.VolumeElement (i); - - int nelfaces = GetNFaces (el.GetType()); - const ELEMENT_FACE * elfaces = GetFaces (el.GetType()); - - for (j = 0; j < 6; j++) - faces.Elem(i)[j] = 0; - for (j = 0; j < nelfaces; j++) - if (elfaces[j][3] == 0) - - { // triangle - - int facenum; - int facedir; - - INDEX_3 face(el.PNum(elfaces[j][0]), - el.PNum(elfaces[j][1]), - el.PNum(elfaces[j][2])); - - facedir = 0; - if (face.I1() > face.I2()) - { - swap (face.I1(), face.I2()); - facedir += 1; - } - if (face.I2() > face.I3()) - { - swap (face.I2(), face.I3()); - facedir += 2; - } - if (face.I1() > face.I2()) - { - swap (face.I1(), face.I2()); - facedir += 4; - } - - if (vert2face.Used (face)) - facenum = vert2face.Get(face); - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); - face2vert.Append (hface); - // face2vert.SetSize(face2vert.Size()+1); - } - - faces.Elem(i)[j] = 8*(facenum-1)+facedir+1; - } - - else - - { - // quad - int facenum; - int facedir; - INDEX_4Q face4(el.PNum(elfaces[j][0]), - el.PNum(elfaces[j][1]), - el.PNum(elfaces[j][2]), - el.PNum(elfaces[j][3])); - - facedir = 0; - if (min2 (face4.I1(), face4.I2()) > - min2 (face4.I4(), face4.I3())) - { // z - flip - facedir += 1; - swap (face4.I1(), face4.I4()); - swap (face4.I2(), face4.I3()); - } - if (min2 (face4.I1(), face4.I4()) > - min2 (face4.I2(), face4.I3())) - { // x - flip - facedir += 2; - swap (face4.I1(), face4.I2()); - swap (face4.I3(), face4.I4()); - } - if (face4.I2() > face4.I4()) - { // diagonal flip - facedir += 4; - swap (face4.I2(), face4.I4()); - } - // face4.Sort(); - - INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); - - if (vert2face.Used (face)) - { - facenum = vert2face.Get(face); - } - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - // face2vert.SetSize(face2vert.Size()+1); - - INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I4()); - face2vert.Append (hface); - } - - faces.Elem(i)[j] = 8*(facenum-1)+facedir+1; - } - } - - face2surfel.SetSize(nfa+nse); - for (i = 1; i <= face2surfel.Size(); i++) - face2surfel.Elem(i) = 0; - - for (i = 1; i <= nse; i++) - { - const Element2d & el = mesh.SurfaceElement (i); - - const ELEMENT_FACE * elfaces = GetFaces (el.GetType()); - - if (elfaces[0][3] == 0) - - { // triangle - - int facenum; - int facedir; - - INDEX_3 face(el.PNum(elfaces[0][0]), - el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2])); - - facedir = 0; - if (face.I1() > face.I2()) - { - swap (face.I1(), face.I2()); - facedir += 1; - } - if (face.I2() > face.I3()) - { - swap (face.I2(), face.I3()); - facedir += 2; - } - if (face.I1() > face.I2()) - { - swap (face.I1(), face.I2()); - facedir += 4; - } - - if (vert2face.Used (face)) - facenum = vert2face.Get(face); - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - // face2vert.SetSize(face2vert.Size()+1); - INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); - face2vert.Append (hface); - } - - surffaces.Elem(i) = 8*(facenum-1)+facedir+1; - face2surfel.Elem(facenum) = i; - } - - else - - { - // quad - int facenum; - int facedir; - - INDEX_4Q face4(el.PNum(elfaces[0][0]), - el.PNum(elfaces[0][1]), - el.PNum(elfaces[0][2]), - el.PNum(elfaces[0][3])); - - facedir = 0; - if (min2 (face4.I1(), face4.I2()) > - min2 (face4.I4(), face4.I3())) - { // z - orientation - facedir += 1; - swap (face4.I1(), face4.I4()); - swap (face4.I2(), face4.I3()); - } - if (min2 (face4.I1(), face4.I4()) > - min2 (face4.I2(), face4.I3())) - { // x - orientation - facedir += 2; - swap (face4.I1(), face4.I2()); - swap (face4.I3(), face4.I4()); - } - if (face4.I2() > face4.I4()) - { - facedir += 4; - swap (face4.I2(), face4.I4()); - } - - INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); - - if (vert2face.Used (face)) - facenum = vert2face.Get(face); - else - { - nfa++; - vert2face.Set (face, nfa); - facenum = nfa; - - // face2vert.SetSize(face2vert.Size()+1); - INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3()); - face2vert.Append (hface); - /* - face2vert.Last()[0] = face4.I1(); - face2vert.Last()[1] = face4.I2(); - face2vert.Last()[2] = face4.I3(); - face2vert.Last()[3] = face4.I3(); - */ - } - - surffaces.Elem(i) = 8*(facenum-1)+facedir+1; - face2surfel.Elem(facenum) = i; - } - } - - - surf2volelement.SetSize (nse); - for (i = 1; i <= nse; i++) - { - surf2volelement.Elem(i)[0] = 0; - surf2volelement.Elem(i)[1] = 0; - } - for (i = 1; i <= ne; i++) - for (j = 0; j < 6; j++) + // ensure all coarse grid and intermediate level edges + cnt = 0; + for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) { - int fnum = (faces.Get(i)[j]+7) / 8; - if (fnum > 0 && face2surfel.Elem(fnum)) + /* JS, Oct 2009 + int pa[2]; + pa[0] = mesh.mlbetweennodes[i].I1(); + pa[1] = mesh.mlbetweennodes[i].I2(); + if (pa[0] > pa[1]) Swap (pa[0], pa[1]); + if (pa[0] > 0) + cnt.Elem(pa[0])++; + */ + INDEX_2 parents = mesh.mlbetweennodes[i]; + parents.Sort(); + if (parents[0] >= PointIndex::BASE) cnt[parents[0]]++; + } + + TABLE vert2vertcoarse (cnt); + for (int i = mesh.mlbetweennodes.Begin(); i < mesh.mlbetweennodes.End(); i++) + { + /* + int pa[2]; + pa[0] = mesh.mlbetweennodes[i].I1(); + pa[1] = mesh.mlbetweennodes[i].I2(); + if (pa[0] > pa[1]) swap (pa[0], pa[1]); + if (pa[0] > 0) + vert2vertcoarse.AddSave1 (pa[0], pa[1]); + */ + INDEX_2 parents = mesh.mlbetweennodes[i]; + parents.Sort(); + if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]); + } + + + Array edgenr(nv); + Array edgeflag(nv); + edgeflag = 0; + + ned = edge2vert.Size(); + Array missing; + + for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) + { + for (int j = 0; j < vert2edge[i].Size(); j++) { - int sel = face2surfel.Elem(fnum); - surf2volelement.Elem(sel)[1] = - surf2volelement.Elem(sel)[0]; - surf2volelement.Elem(sel)[0] = i; + int ednr = vert2edge[i][j]; + int i2 = edge2vert.Get(ednr)[1]; + edgeflag[i2] = i; + edgenr[i2] = ednr; + } + for (int j = 0; j << vert2vertcoarse[i].Size(); j++) + { + int v2 = vert2vertcoarse[i][j]; + if (edgeflag[v2] < i) + { + ned++; + edgenr[v2] = ned; + edgeflag[v2] = i; + missing.Append (INDEX_3(i,v2,ned)); + } + } + + for (int j = 0; j < (*vert2element)[i].Size(); j++) + { + int elnr = (*vert2element)[i][j]; + const Element & el = mesh.VolumeElement (elnr); + + int neledges = GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + if (edge.I1() != i) + continue; + + if (edgeflag[edge.I2()] < i) + { + ned++; + edgenr[edge.I2()] = ned; + edgeflag[edge.I2()] = i; + } + + int edgenum = edgenr[edge.I2()]; + if (edgedir) edgenum *= -1; + edges.Elem(elnr)[k] = edgenum; + } + } + + for (int j = 0; j < (*vert2surfelement)[i].Size(); j++) + { + int elnr = (*vert2surfelement)[i][j]; + const Element2d & el = mesh.SurfaceElement (elnr); + + int neledges = GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + if (edge.I1() != i) + continue; + + if (edgeflag[edge.I2()] < i) + { + ned++; + edgenr[edge.I2()] = ned; + edgeflag[edge.I2()] = i; + } + + int edgenum = edgenr[edge.I2()]; + if (edgedir) edgenum *= -1; + surfedges.Elem(elnr)[k] = edgenum; + } + } + + for (int j = 0; j < (*vert2segment)[i].Size(); j++) + { + int elnr = (*vert2segment)[i][j]; + const Segment & el = mesh.LineSegment (elnr); + + INDEX_2 edge(el[0], el[1]); + + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + if (edge.I1() != i) + continue; + + if (edgeflag[edge.I2()] < i) + { + ned++; + edgenr[edge.I2()] = ned; + edgeflag[edge.I2()] = i; + } + int edgenum = edgenr[edge.I2()]; + + if (edgedir) edgenum *= -1; + segedges.Elem(elnr) = edgenum; } } - face2vert.SetAllocSize (face2vert.Size()); - /* - *testout << "face2vert: "; - face2vert.PrintMemInfo(cout); - *testout << "faces: "; - faces.PrintMemInfo(cout); - *testout << "hashtable: "; - vert2face.PrintMemInfo(cout); - */ + edge2vert.SetSize (ned); + for (int i = 1; i <= ne; i++) + { + const Element & el = mesh.VolumeElement (i); + + int neledges = GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); -#ifdef PARALLEL - (*testout) << " RESET Paralleltop" << endl; + int edgenum = abs (edges.Elem(i)[k]); - paralleltop.Reset (); -#endif + edge2vert.Elem(edgenum)[0] = edge.I1(); + edge2vert.Elem(edgenum)[1] = edge.I2(); + } + } + for (int i = 1; i <= nse; i++) + { + const Element2d & el = mesh.SurfaceElement (i); + + int neledges = GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); - Array face_els(nfa), face_surfels(nfa); - face_els = 0; - face_surfels = 0; - Array hfaces; - for (i = 1; i <= ne; i++) + int edgenum = abs (surfedges.Elem(i)[k]); + + edge2vert.Elem(edgenum)[0] = edge.I1(); + edge2vert.Elem(edgenum)[1] = edge.I2(); + } + } + + for (int i = 1; i <= nseg; i++) + { + const Segment & el = mesh.LineSegment (i); + + INDEX_2 edge(el[0], el[1]); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + int edgenum = abs (segedges.Elem(i)); + + edge2vert.Elem(edgenum)[0] = edge.I1(); + edge2vert.Elem(edgenum)[1] = edge.I2(); + } + + for (int i = 1; i <= missing.Size(); i++) + { + INDEX_3 i3 = missing.Get(i); + edge2vert.Elem(i3.I3())[0] = i3.I1(); + edge2vert.Elem(i3.I3())[1] = i3.I2(); + } + + + /* + (*testout) << "edge table:" << endl; + (*testout) << "edge2vert:" << endl; + for (int i = 1; i <= edge2vert.Size(); i++) + (*testout) << "edge " << i << ", v1,2 = " << edge2vert.Elem(i)[0] << ", " << edge2vert.Elem(i)[1] << endl; + (*testout) << "surfedges:" << endl; + for (int i = 1; i <= surfedges.Size(); i++) + (*testout) << "el " << i << ", edges = " + << surfedges.Elem(i)[0] << ", " + << surfedges.Elem(i)[1] << ", " + << surfedges.Elem(i)[2] << endl; + */ + + + } + + + // cout << "build edges done" << endl; + + // generate faces + if (buildfaces) // && mesh.GetDimension() == 3) + { + static int timer2 = NgProfiler::CreateTimer ("topology::buildfaces"); + NgProfiler::RegionTimer reg2 (timer2); + + PrintMessage (5, "Update faces "); + + + + + + /* + + + faces.SetSize(ne); + surffaces.SetSize(nse); + + // face2vert.SetSize(0); // keep old faces + nfa = face2vert.Size(); + // INDEX_3_HASHTABLE vert2face(ne+nse+1); + INDEX_3_CLOSED_HASHTABLE vert2face(8*ne+2*nse+nfa+2); + + for (int i = 1; i <= face2vert.Size(); i++) { - GetElementFaces (i, hfaces); - for (j = 0; j < hfaces.Size(); j++) - face_els[hfaces[j]-1]++; + INDEX_3 f; + f.I1() = face2vert.Get(i)[0]; + f.I2() = face2vert.Get(i)[1]; + f.I3() = face2vert.Get(i)[2]; + vert2face.Set (f, i); } - for (i = 1; i <= nse; i++) - face_surfels[GetSurfaceElementFace (i)-1]++; - - if (ne) + + for (int i = 1; i <= ne; i++) { - int cnt_err = 0; - for (i = 0; i < nfa; i++) - { - /* - (*testout) << "face " << i << " has " << int(face_els[i]) << " els, " - << int(face_surfels[i]) << " surfels, tot = " - << face_els[i] + face_surfels[i] << endl; - */ - if (face_els[i] + face_surfels[i] == 1) - { - cnt_err++; -#ifdef PARALLEL - if ( ntasks > 1 ) - { - if ( !paralleltop.DoCoarseUpdate() ) continue; - // "illegal" faces are exchange faces - /* - (*testout) << "exchange face : " << i << endl; - (*testout) << "points = " << face2vert[i] << endl; - */ -// (*testout) << "global points = "; -// for ( int j = 0; j < 3; j++ ) -// // (*testout) << face2vert[i].I(j+1) << " -> " -// // << paralleltop.GetLoc2Glob_Vert( face2vert[i].I(j+1) ) << ", "; -// (*testout) << endl; -// if ( !paralleltop.IsExchangeFace (i+1) ) -// paralleltop.SetRefinementFace (i+1); + const Element & el = mesh.VolumeElement (i); + + int nelfaces = GetNFaces (el.GetType()); + const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); + + for (int j = 0; j < 6; j++) + faces.Elem(i)[j] = 0; + for (int j = 0; j < nelfaces; j++) + if (elfaces[j][3] == 0) + + { // triangle + + int facenum; + int facedir; + + INDEX_3 face(el.PNum(elfaces[j][0]), + el.PNum(elfaces[j][1]), + el.PNum(elfaces[j][2])); + + facedir = 0; + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 1; + } + if (face.I2() > face.I3()) + { + swap (face.I2(), face.I3()); + facedir += 2; + } + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 4; + } + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); + face2vert.Append (hface); + // face2vert.SetSize(face2vert.Size()+1); + } + + faces.Elem(i)[j] = 8*(facenum-1)+facedir+1; + } + + else + + { + // quad + int facenum; + int facedir; + INDEX_4Q face4(el.PNum(elfaces[j][0]), + el.PNum(elfaces[j][1]), + el.PNum(elfaces[j][2]), + el.PNum(elfaces[j][3])); + + facedir = 0; + if (min2 (face4.I1(), face4.I2()) > + min2 (face4.I4(), face4.I3())) + { // z - flip + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - flip + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { // diagonal flip + facedir += 4; + swap (face4.I2(), face4.I4()); + } + // face4.Sort(); + + INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); + + if (vert2face.Used (face)) + { + facenum = vert2face.Get(face); + } + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; - paralleltop.SetExchangeFace (i+1); - - for (int j = 0; j < 4; j++) - { - if ( face2vert[i].I(j+1) > 0 ) - paralleltop.SetExchangeVert(face2vert[i].I(j+1)); - } - - Array faceedges; - GetFaceEdges (i+1, faceedges); - for ( int j = 0; j < faceedges.Size(); j++) - { - paralleltop.SetExchangeEdge ( faceedges[j] ); - int v1, v2; - GetEdgeVertices(faceedges[j], v1, v2 ); - } - - /* - (*testout) << "pos = "; - for (int j = 0; j < 4; j++) - if (face2vert[i].I(j+1) >= 1) - (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " "; - (*testout) << endl; - */ - } - else - { -#endif - (*testout) << "illegal face : " << i << endl; - (*testout) << "points = " << face2vert[i] << endl; - (*testout) << "pos = "; - for (j = 0; j < 4; j++) - if (face2vert[i].I(j+1) >= 1) - (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " "; - (*testout) << endl; + // face2vert.SetSize(face2vert.Size()+1); - FlatArray vertels = GetVertexElements (face2vert[i].I(1)); - for (int k = 0; k < vertels.Size(); k++) - { - int elfaces[10], orient[10]; - int nf = GetElementFaces (vertels[k], elfaces, orient); - for (int l = 0; l < nf; l++) - if (elfaces[l] == i) + INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I4()); + face2vert.Append (hface); + } + + faces.Elem(i)[j] = 8*(facenum-1)+facedir+1; + } + } + + face2surfel.SetSize(nfa+nse); + for (int i = 1; i <= face2surfel.Size(); i++) + face2surfel.Elem(i) = 0; + + for (int i = 1; i <= nse; i++) + { + const Element2d & el = mesh.SurfaceElement (i); + + const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); + + if (elfaces[0][3] == 0) + + { // triangle + + int facenum; + int facedir; + + INDEX_3 face(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2])); + + facedir = 0; + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 1; + } + if (face.I2() > face.I3()) + { + swap (face.I2(), face.I3()); + facedir += 2; + } + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 4; + } + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + // face2vert.SetSize(face2vert.Size()+1); + INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); + face2vert.Append (hface); + } + + surffaces.Elem(i) = 8*(facenum-1)+facedir+1; + face2surfel.Elem(facenum) = i; + } + + else + + { + // quad + int facenum; + int facedir; + + INDEX_4Q face4(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2]), + el.PNum(elfaces[0][3])); + + facedir = 0; + if (min2 (face4.I1(), face4.I2()) > + min2 (face4.I4(), face4.I3())) + { // z - orientation + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - orientation + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + // face2vert.SetSize(face2vert.Size()+1); + INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3()); + face2vert.Append (hface); + } + + surffaces.Elem(i) = 8*(facenum-1)+facedir+1; + face2surfel.Elem(facenum) = i; + } + } + + + surf2volelement.SetSize (nse); + for (int i = 1; i <= nse; i++) + { + surf2volelement.Elem(i)[0] = 0; + surf2volelement.Elem(i)[1] = 0; + } + for (int i = 1; i <= ne; i++) + for (int j = 0; j < 6; j++) + { + int fnum = (faces.Get(i)[j]+7) / 8; + if (fnum > 0 && face2surfel.Elem(fnum)) + { + int sel = face2surfel.Elem(fnum); + surf2volelement.Elem(sel)[1] = + surf2volelement.Elem(sel)[0]; + surf2volelement.Elem(sel)[0] = i; + } + } + + face2vert.SetAllocSize (face2vert.Size()); + */ + + + + + + + + faces.SetSize(ne); + surffaces.SetSize(nse); + + int oldnfa = face2vert.Size(); + + cnt = 0; + for (int i = 0; i < face2vert.Size(); i++) + cnt[face2vert[i][0]]++; + TABLE vert2oldface(cnt); + for (int i = 0; i < face2vert.Size(); i++) + vert2oldface.AddSave (face2vert[i][0], i); + + + for (int elnr = 0; elnr < ne; elnr++) + for (int j = 0; j < 6; j++) + faces[elnr][j] = 0; + + + int max_face_on_vertex = 0; + for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) + { + int onv = vert2oldface[i].Size() + (*vert2element)[i].Size() + (*vert2surfelement)[i].Size(); + max_face_on_vertex = max (onv, max_face_on_vertex); + } + + + for (int pass = 1; pass <= 2; pass++) + { + nfa = oldnfa; + for (int v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) + { + INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); + + for (int j = 0; j < vert2oldface[v].Size(); j++) + { + int fnr = vert2oldface[v][j]; + INDEX_3 face (face2vert[fnr].I1(), + face2vert[fnr].I2(), + face2vert[fnr].I3()); + vert2face.Set (face, fnr+1); + } + + // cout << "inherited faces: " << endl << vert2face << endl; + + + for (int j = 0; j < (*vert2element)[v].Size(); j++) + { + int elnr = (*vert2element)[v][j]; + const Element & el = mesh.VolumeElement (elnr); + + int nelfaces = GetNFaces (el.GetType()); + const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); + + for (int j = 0; j < nelfaces; j++) + if (elfaces[j][3] == 0) + + { // triangle + + int facenum; + int facedir; + + INDEX_3 face(el.PNum(elfaces[j][0]), + el.PNum(elfaces[j][1]), + el.PNum(elfaces[j][2])); + + facedir = 0; + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 1; + } + if (face.I2() > face.I3()) + { + swap (face.I2(), face.I3()); + facedir += 2; + } + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 4; + } + + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); + if (pass == 2) face2vert.Append (hface); + } + // cout << "face = " << face << " elnr = " << elnr << endl; + faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1; + } + + else + + { + // quad + int facenum; + int facedir; + INDEX_4Q face4(el.PNum(elfaces[j][0]), + el.PNum(elfaces[j][1]), + el.PNum(elfaces[j][2]), + el.PNum(elfaces[j][3])); + + facedir = 0; + if (min2 (face4.I1(), face4.I2()) > + min2 (face4.I4(), face4.I3())) + { // z - flip + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - flip + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { // diagonal flip + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + + INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); + + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + { + facenum = vert2face.Get(face); + } + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I4()); + if (pass == 2) face2vert.Append (hface); + } + + faces.Elem(elnr)[j] = 8*(facenum-1)+facedir+1; + } + } + + for (int j = 0; j < (*vert2surfelement)[v].Size(); j++) + { + int elnr = (*vert2surfelement)[v][j]; + // cout << "surfelnr = " << elnr << endl; + const Element2d & el = mesh.SurfaceElement (elnr); + + const ELEMENT_FACE * elfaces = GetFaces1 (el.GetType()); + + if (elfaces[0][3] == 0) + + { // triangle + + int facenum; + int facedir; + + INDEX_3 face(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2])); + + // cout << "face = " << face << endl; + + facedir = 0; + if (face.I1() > face.I2()) { - (*testout) << "is face of element " << vertels[k] << endl; - - if (mesh.coarsemesh && mesh.hpelements->Size() == mesh.GetNE() ) - { - const HPRefElement & hpref_el = - (*mesh.hpelements) [ mesh.VolumeElement (vertels[k]).hp_elnr]; - (*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl; - } - + swap (face.I1(), face.I2()); + facedir += 1; } - } -#ifdef PARALLEL - } -#endif + if (face.I2() > face.I3()) + { + swap (face.I2(), face.I3()); + facedir += 2; + } + if (face.I1() > face.I2()) + { + swap (face.I1(), face.I2()); + facedir += 4; + } + + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face.I1(),face.I2(),face.I3(),0); + if (pass == 2) face2vert.Append (hface); + } + + // cout << "face = " << face << " selnr = " << elnr << endl; + surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1; + // face2surfel.Elem(facenum) = elnr; + } + + else + + { + // quad + int facenum; + int facedir; + + INDEX_4Q face4(el.PNum(elfaces[0][0]), + el.PNum(elfaces[0][1]), + el.PNum(elfaces[0][2]), + el.PNum(elfaces[0][3])); + + facedir = 0; + if (min2 (face4.I1(), face4.I2()) > + min2 (face4.I4(), face4.I3())) + { // z - orientation + facedir += 1; + swap (face4.I1(), face4.I4()); + swap (face4.I2(), face4.I3()); + } + if (min2 (face4.I1(), face4.I4()) > + min2 (face4.I2(), face4.I3())) + { // x - orientation + facedir += 2; + swap (face4.I1(), face4.I2()); + swap (face4.I3(), face4.I4()); + } + if (face4.I2() > face4.I4()) + { + facedir += 4; + swap (face4.I2(), face4.I4()); + } + + INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); + if (face.I1() != v) continue; + + if (vert2face.Used (face)) + facenum = vert2face.Get(face); + else + { + nfa++; + vert2face.Set (face, nfa); + facenum = nfa; + + INDEX_4 hface(face4.I1(),face4.I2(),face4.I3(),face4.I3()); + if (pass == 2) face2vert.Append (hface); + } + + surffaces.Elem(elnr) = 8*(facenum-1)+facedir+1; + // face2surfel.Elem(facenum) = elnr; + } + } + } + face2vert.SetAllocSize (nfa); + } + + + + face2surfel.SetSize (nfa); + face2surfel = 0; + for (int i = 1; i <= nse; i++) + face2surfel.Elem(GetSurfaceElementFace(i)) = i; + + /* + cout << "build table complete" << endl; + + cout << "faces = " << endl; + + cout << "face2vert = " << endl << face2vert << endl; + cout << "surffaces = " << endl << surffaces << endl; + cout << "face2surfel = " << endl << face2surfel << endl; + */ + + + surf2volelement.SetSize (nse); + for (int i = 1; i <= nse; i++) + { + surf2volelement.Elem(i)[0] = 0; + surf2volelement.Elem(i)[1] = 0; + } + for (int i = 1; i <= ne; i++) + for (int j = 0; j < 6; j++) + { + int fnum = (faces.Get(i)[j]+7) / 8; + if (fnum > 0 && face2surfel.Elem(fnum)) + { + int sel = face2surfel.Elem(fnum); + surf2volelement.Elem(sel)[1] = + surf2volelement.Elem(sel)[0]; + surf2volelement.Elem(sel)[0] = i; } } -#ifndef PARALLEL - if (cnt_err) - cout << cnt_err << " elements are not matching !!!" << endl; -#else - if (cnt_err && ntasks == 1) - cout << cnt_err << " elements are not matching !!!" << endl; - else if (cnt_err && ntasks > 1) - { - cout << "p" << id << ": " << cnt_err << " elements are not local" << endl; - isparallel = 1; - } - else if ( ntasks > 1 ) - cout << "p" << id << ": " << "Partition " << id << " is totally local" << endl; + face2vert.SetAllocSize (face2vert.Size()); + + + + + + + + + + + + + + + + + + // face table complete + + + + /* + *testout << "face2vert: "; + face2vert.PrintMemInfo(cout); + *testout << "faces: "; + faces.PrintMemInfo(cout); + *testout << "hashtable: "; + vert2face.PrintMemInfo(cout); + */ + +#ifdef PARALLEL + (*testout) << " RESET Paralleltop" << endl; + + paralleltop.Reset (); #endif - } + Array face_els(nfa), face_surfels(nfa); + face_els = 0; + face_surfels = 0; + Array hfaces; + for (int i = 1; i <= ne; i++) + { + GetElementFaces (i, hfaces); + for (int j = 0; j < hfaces.Size(); j++) + face_els[hfaces[j]-1]++; + } + for (int i = 1; i <= nse; i++) + face_surfels[GetSurfaceElementFace (i)-1]++; + + + if (ne) + { + int cnt_err = 0; + for (int i = 0; i < nfa; i++) + { + /* + (*testout) << "face " << i << " has " << int(face_els[i]) << " els, " + << int(face_surfels[i]) << " surfels, tot = " + << face_els[i] + face_surfels[i] << endl; + */ + if (face_els[i] + face_surfels[i] == 1) + { + cnt_err++; +#ifdef PARALLEL + if ( ntasks > 1 ) + { + if ( !paralleltop.DoCoarseUpdate() ) continue; + // "illegal" faces are exchange faces + /* + (*testout) << "exchange face : " << i << endl; + (*testout) << "points = " << face2vert[i] << endl; + */ + // (*testout) << "global points = "; + // for ( int j = 0; j < 3; j++ ) + // // (*testout) << face2vert[i].I(j+1) << " -> " + // // << paralleltop.GetLoc2Glob_Vert( face2vert[i].I(j+1) ) << ", "; + // (*testout) << endl; + // if ( !paralleltop.IsExchangeFace (i+1) ) + // paralleltop.SetRefinementFace (i+1); + + paralleltop.SetExchangeFace (i+1); + + for (int j = 0; j < 4; j++) + { + if ( face2vert[i].I(j+1) > 0 ) + paralleltop.SetExchangeVert(face2vert[i].I(j+1)); + } + + Array faceedges; + GetFaceEdges (i+1, faceedges); + for ( int j = 0; j < faceedges.Size(); j++) + { + paralleltop.SetExchangeEdge ( faceedges[j] ); + int v1, v2; + GetEdgeVertices(faceedges[j], v1, v2 ); + } + + /* + (*testout) << "pos = "; + for (int j = 0; j < 4; j++) + if (face2vert[i].I(j+1) >= 1) + (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " "; + (*testout) << endl; + */ + } + else + { +#endif + (*testout) << "illegal face : " << i << endl; + (*testout) << "points = " << face2vert[i] << endl; + (*testout) << "pos = "; + for (int j = 0; j < 4; j++) + if (face2vert[i].I(j+1) >= 1) + (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " "; + (*testout) << endl; + + FlatArray vertels = GetVertexElements (face2vert[i].I(1)); + for (int k = 0; k < vertels.Size(); k++) + { + int elfaces[10], orient[10]; + int nf = GetElementFaces (vertels[k], elfaces, orient); + for (int l = 0; l < nf; l++) + if (elfaces[l] == i) + { + (*testout) << "is face of element " << vertels[k] << endl; + + if (mesh.coarsemesh && mesh.hpelements->Size() == mesh.GetNE() ) + { + const HPRefElement & hpref_el = + (*mesh.hpelements) [ mesh.VolumeElement (vertels[k]).hp_elnr]; + (*testout) << "coarse eleme = " << hpref_el.coarse_elnr << endl; + } + + } + } +#ifdef PARALLEL + } +#endif + } + } + +#ifndef PARALLEL + if (cnt_err) + cout << cnt_err << " elements are not matching !!!" << endl; +#else + if (cnt_err && ntasks == 1) + cout << cnt_err << " elements are not matching !!!" << endl; + else if (cnt_err && ntasks > 1) + { + cout << "p" << id << ": " << cnt_err << " elements are not local" << endl; + isparallel = 1; + } + else if ( ntasks > 1 ) + cout << "p" << id << ": " << "Partition " << id << " is totally local" << endl; +#endif + + } #ifdef PARALLEL - if ( isparallel ) - { - paralleltop.Update(); - if ( paralleltop.DoCoarseUpdate() ) - { - paralleltop.UpdateCoarseGrid(); - } - else - { - // paralleltop.UpdateRefinement(); - } - // paralleltop.Print(); - } + if ( isparallel ) + { + paralleltop.Update(); + if ( paralleltop.DoCoarseUpdate() ) + { + paralleltop.UpdateCoarseGrid(); + } + else + { + // paralleltop.UpdateRefinement(); + } + // paralleltop.Print(); + } #endif @@ -815,127 +1158,127 @@ void MeshTopology :: Update() - } + } - /* -for (i = 1; i <= ne; i++) - { - (*testout) << "Element " << i << endl; - (*testout) << "PNums " << endl; - for( int l=1;l<=8;l++) *testout << mesh.VolumeElement(i).PNum(l) << "\t"; - *testout << endl; - (*testout) << "edges: " << endl; - for (j = 0; j < 9; j++) - (*testout) << edges.Elem(i)[j] << " "; - (*testout) << "faces: " << endl; - for (j = 0; j < 6; j++)m - (*testout) << faces.Elem(i)[j] << " "; - } + /* + for (i = 1; i <= ne; i++) + { + (*testout) << "Element " << i << endl; + (*testout) << "PNums " << endl; + for( int l=1;l<=8;l++) *testout << mesh.VolumeElement(i).PNum(l) << "\t"; + *testout << endl; + (*testout) << "edges: " << endl; + for (j = 0; j < 9; j++) + (*testout) << edges.Elem(i)[j] << " "; + (*testout) << "faces: " << endl; + for (j = 0; j < 6; j++)m + (*testout) << faces.Elem(i)[j] << " "; + } - for (i = 1; i <= nse; i++) - { - (*testout) << "SElement " << i << endl; - (*testout) << "PNums " << endl; - for( int l=1;l<=4;l++) *testout << mesh.SurfaceElement(i).PNum(l) << "\t"; - *testout << endl; - } - */ - timestamp = NextTimeStamp(); -} + for (i = 1; i <= nse; i++) + { + (*testout) << "SElement " << i << endl; + (*testout) << "PNums " << endl; + for( int l=1;l<=4;l++) *testout << mesh.SurfaceElement(i).PNum(l) << "\t"; + *testout << endl; + } + */ + timestamp = NextTimeStamp(); + } -const Point3d * MeshTopology :: GetVertices (ELEMENT_TYPE et) -{ - static Point3d segm_points [] = - { Point3d (1, 0, 0), - Point3d (0, 0, 0) }; + const Point3d * MeshTopology :: GetVertices (ELEMENT_TYPE et) + { + static Point3d segm_points [] = + { Point3d (1, 0, 0), + Point3d (0, 0, 0) }; - static Point3d trig_points [] = - { Point3d ( 1, 0, 0 ), - Point3d ( 0, 1, 0 ), - Point3d ( 0, 0, 0 ) }; + static Point3d trig_points [] = + { Point3d ( 1, 0, 0 ), + Point3d ( 0, 1, 0 ), + Point3d ( 0, 0, 0 ) }; - static Point3d quad_points [] = - { Point3d ( 0, 0, 0 ), - Point3d ( 1, 0, 0 ), - Point3d ( 1, 1, 0 ), - Point3d ( 0, 1, 0 ) }; + static Point3d quad_points [] = + { Point3d ( 0, 0, 0 ), + Point3d ( 1, 0, 0 ), + Point3d ( 1, 1, 0 ), + Point3d ( 0, 1, 0 ) }; - static Point3d tet_points [] = - { Point3d ( 1, 0, 0 ), - Point3d ( 0, 1, 0 ), - Point3d ( 0, 0, 1 ), - Point3d ( 0, 0, 0 ) }; + static Point3d tet_points [] = + { Point3d ( 1, 0, 0 ), + Point3d ( 0, 1, 0 ), + Point3d ( 0, 0, 1 ), + Point3d ( 0, 0, 0 ) }; - static Point3d pyramid_points [] = - { - Point3d ( 0, 0, 0 ), - Point3d ( 1, 0, 0 ), - Point3d ( 1, 1, 0 ), - Point3d ( 0, 1, 0 ), - Point3d ( 0, 0, 1-1e-7 ), - }; + static Point3d pyramid_points [] = + { + Point3d ( 0, 0, 0 ), + Point3d ( 1, 0, 0 ), + Point3d ( 1, 1, 0 ), + Point3d ( 0, 1, 0 ), + Point3d ( 0, 0, 1-1e-7 ), + }; - static Point3d prism_points[] = - { - Point3d ( 1, 0, 0 ), - Point3d ( 0, 1, 0 ), - Point3d ( 0, 0, 0 ), - Point3d ( 1, 0, 1 ), - Point3d ( 0, 1, 1 ), - Point3d ( 0, 0, 1 ) - }; + static Point3d prism_points[] = + { + Point3d ( 1, 0, 0 ), + Point3d ( 0, 1, 0 ), + Point3d ( 0, 0, 0 ), + Point3d ( 1, 0, 1 ), + Point3d ( 0, 1, 1 ), + Point3d ( 0, 0, 1 ) + }; - static Point3d hex_points [] = - { Point3d ( 0, 0, 0 ), - Point3d ( 1, 0, 0 ), - Point3d ( 1, 1, 0 ), - Point3d ( 0, 1, 0 ), - Point3d ( 0, 0, 1 ), - Point3d ( 1, 0, 1 ), - Point3d ( 1, 1, 1 ), - Point3d ( 0, 1, 1 ) }; + static Point3d hex_points [] = + { Point3d ( 0, 0, 0 ), + Point3d ( 1, 0, 0 ), + Point3d ( 1, 1, 0 ), + Point3d ( 0, 1, 0 ), + Point3d ( 0, 0, 1 ), + Point3d ( 1, 0, 1 ), + Point3d ( 1, 1, 1 ), + Point3d ( 0, 1, 1 ) }; - switch (et) - { - case SEGMENT: - case SEGMENT3: - return segm_points; + switch (et) + { + case SEGMENT: + case SEGMENT3: + return segm_points; - case TRIG: - case TRIG6: - return trig_points; + case TRIG: + case TRIG6: + return trig_points; - case QUAD: - case QUAD6: - case QUAD8: - return quad_points; + case QUAD: + case QUAD6: + case QUAD8: + return quad_points; - case TET: - case TET10: - return tet_points; + case TET: + case TET10: + return tet_points; - case PYRAMID: - return pyramid_points; + case PYRAMID: + return pyramid_points; - case PRISM: - case PRISM12: - return prism_points; + case PRISM: + case PRISM12: + return prism_points; - case HEX: - return hex_points; - default: - cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; - } - return 0; -} + case HEX: + return hex_points; + default: + cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; + } + return 0; + } @@ -944,404 +1287,406 @@ const Point3d * MeshTopology :: GetVertices (ELEMENT_TYPE et) -void MeshTopology :: GetElementEdges (int elnr, Array & eledges) const -{ - int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); - eledges.SetSize (ned); - for (int i = 0; i < ned; i++) - eledges[i] = abs (edges.Get(elnr)[i]); -} -void MeshTopology :: GetElementFaces (int elnr, Array & elfaces, bool withorientation) const -{ - int i; - int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); - elfaces.SetSize (nfa); - for (i = 1; i <= nfa; i++) - { - elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; - if(withorientation) + void MeshTopology :: GetElementEdges (int elnr, Array & eledges) const + { + int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); + eledges.SetSize (ned); + for (int i = 0; i < ned; i++) + eledges[i] = abs (edges.Get(elnr)[i]); + } + void MeshTopology :: GetElementFaces (int elnr, Array & elfaces, bool withorientation) const + { + int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); + elfaces.SetSize (nfa); + + if (!withorientation) + + for (int i = 1; i <= nfa; i++) { + elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; + } + + else + + for (int i = 1; i <= nfa; i++) + { + elfaces.Elem(i) = (faces.Get(elnr)[i-1]-1) / 8 + 1; int orient = (faces.Get(elnr)[i-1]-1) % 8; if(orient == 1 || orient == 2 || orient == 4 || orient == 7) elfaces.Elem(i) *= -1; } - } -} + } -void MeshTopology :: GetElementEdgeOrientations (int elnr, Array & eorient) const -{ - int i; - int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); - eorient.SetSize (ned); - for (i = 1; i <= ned; i++) - eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; -} + void MeshTopology :: GetElementEdgeOrientations (int elnr, Array & eorient) const + { + int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); + eorient.SetSize (ned); + for (int i = 1; i <= ned; i++) + eorient.Elem(i) = (edges.Get(elnr)[i-1] > 0) ? 1 : -1; + } -void MeshTopology :: GetElementFaceOrientations (int elnr, Array & forient) const -{ - int i; - int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); - forient.SetSize (nfa); - for (i = 1; i <= nfa; i++) - forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; -} + void MeshTopology :: GetElementFaceOrientations (int elnr, Array & forient) const + { + int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); + forient.SetSize (nfa); + for (int i = 1; i <= nfa; i++) + forient.Elem(i) = (faces.Get(elnr)[i-1]-1) % 8; + } -int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) const -{ - int i; - // int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); + int MeshTopology :: GetElementEdges (int elnr, int * eledges, int * orient) const + { + // int ned = GetNEdges (mesh.VolumeElement(elnr).GetType()); - if (mesh.GetDimension()==3 || 1) - { - if (orient) - { - for (i = 0; i < 12; i++) - { - if (!edges.Get(elnr)[i]) return i; - eledges[i] = abs (edges.Get(elnr)[i]); - orient[i] = (edges.Get(elnr)[i] > 0 ) ? 1 : -1; - } - } - else - { - for (i = 0; i < 12; i++) - { - if (!edges.Get(elnr)[i]) return i; - eledges[i] = abs (edges.Get(elnr)[i]); - } - } - return 12; - } - else - { - throw NgException("rethink implementation"); - /* - if (orient) - { + if (mesh.GetDimension()==3 || 1) + { + if (orient) + { + for (int i = 0; i < 12; i++) + { + if (!edges.Get(elnr)[i]) return i; + eledges[i] = abs (edges.Get(elnr)[i]); + orient[i] = (edges.Get(elnr)[i] > 0 ) ? 1 : -1; + } + } + else + { + for (int i = 0; i < 12; i++) + { + if (!edges.Get(elnr)[i]) return i; + eledges[i] = abs (edges.Get(elnr)[i]); + } + } + return 12; + } + else + { + throw NgException("rethink implementation"); + /* + if (orient) + { for (i = 0; i < 4; i++) - { - if (!surfedges.Get(elnr)[i]) return i; - eledges[i] = abs (surfedges.Get(elnr)[i]); - orient[i] = (surfedges.Get(elnr)[i] > 0 ) ? 1 : -1; - } - } - else - { + { + if (!surfedges.Get(elnr)[i]) return i; + eledges[i] = abs (surfedges.Get(elnr)[i]); + orient[i] = (surfedges.Get(elnr)[i] > 0 ) ? 1 : -1; + } + } + else + { if (!surfedges.Get(elnr)[i]) return i; for (i = 0; i < 4; i++) - eledges[i] = abs (surfedges.Get(elnr)[i]); - } + eledges[i] = abs (surfedges.Get(elnr)[i]); + } */ - return 4; - // return GetSurfaceElementEdges (elnr, eledges, orient); - } -} + return 4; + // return GetSurfaceElementEdges (elnr, eledges, orient); + } + } -int MeshTopology :: GetElementFaces (int elnr, int * elfaces, int * orient) const -{ - int i; - // int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); - if (orient) - { - for (i = 0; i < 6; i++) - { - if (!faces.Get(elnr)[i]) return i; - elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; - orient[i] = (faces.Get(elnr)[i]-1) % 8; - } - } - else - { - for (i = 0; i < 6; i++) - { - if (!faces.Get(elnr)[i]) return i; - elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; - } - } - return 6; -} + int MeshTopology :: GetElementFaces (int elnr, int * elfaces, int * orient) const + { + // int nfa = GetNFaces (mesh.VolumeElement(elnr).GetType()); + if (orient) + { + for (int i = 0; i < 6; i++) + { + if (!faces.Get(elnr)[i]) return i; + elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; + orient[i] = (faces.Get(elnr)[i]-1) % 8; + } + } + else + { + for (int i = 0; i < 6; i++) + { + if (!faces.Get(elnr)[i]) return i; + elfaces[i] = (faces.Get(elnr)[i]-1) / 8 + 1; + } + } + return 6; + } -void MeshTopology :: GetSurfaceElementEdges (int elnr, Array & eledges) const -{ - int i; - if (mesh.GetDimension()==3 || 1) - { - int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); - eledges.SetSize (ned); - for (i = 1; i <= ned; i++) - eledges.Elem(i) = abs (surfedges.Get(elnr)[i-1]); - } - else - { - cout << "surfeledge(" << elnr << ") = " << flush; - eledges.SetSize(1); - eledges.Elem(1) = abs (segedges.Get(elnr)); - cout << eledges.Elem(1) << endl; - } -} + void MeshTopology :: GetSurfaceElementEdges (int elnr, Array & eledges) const + { + int i; + if (mesh.GetDimension()==3 || 1) + { + int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); + eledges.SetSize (ned); + for (i = 1; i <= ned; i++) + eledges.Elem(i) = abs (surfedges.Get(elnr)[i-1]); + } + else + { + cout << "surfeledge(" << elnr << ") = " << flush; + eledges.SetSize(1); + eledges.Elem(1) = abs (segedges.Get(elnr)); + cout << eledges.Elem(1) << endl; + } + } -int MeshTopology :: GetSurfaceElementFace (int elnr) const -{ - return (surffaces.Get(elnr)-1) / 8 + 1; -} + int MeshTopology :: GetSurfaceElementFace (int elnr) const + { + return (surffaces.Get(elnr)-1) / 8 + 1; + } -void MeshTopology :: -GetSurfaceElementEdgeOrientations (int elnr, Array & eorient) const -{ - int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); - eorient.SetSize (ned); - for (int i = 1; i <= ned; i++) - eorient.Elem(i) = (surfedges.Get(elnr)[i-1] > 0) ? 1 : -1; -} + void MeshTopology :: + GetSurfaceElementEdgeOrientations (int elnr, Array & eorient) const + { + int ned = GetNEdges (mesh.SurfaceElement(elnr).GetType()); + eorient.SetSize (ned); + for (int i = 1; i <= ned; i++) + eorient.Elem(i) = (surfedges.Get(elnr)[i-1] > 0) ? 1 : -1; + } -int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const -{ - return (surffaces.Get(elnr)-1) % 8; -} + int MeshTopology :: GetSurfaceElementFaceOrientation (int elnr) const + { + return (surffaces.Get(elnr)-1) % 8; + } -int MeshTopology :: GetSurfaceElementEdges (int elnr, int * eledges, int * orient) const -{ - int i; - if (mesh.GetDimension() == 3 || 1) - { - if (orient) - { - for (i = 0; i < 4; i++) - { - if (!surfedges.Get(elnr)[i]) return i; - eledges[i] = abs (surfedges.Get(elnr)[i]); - orient[i] = (surfedges.Get(elnr)[i] > 0 ) ? 1 : -1; - } - } - else - { - for (i = 0; i < 4; i++) - { - if (!surfedges.Get(elnr)[i]) return i; - eledges[i] = abs (surfedges.Get(elnr)[i]); - } - } - return 4; - } - else - { - eledges[0] = abs (segedges.Get(elnr)); - if (orient) - orient[0] = segedges.Get(elnr) > 0 ? 1 : -1; - } - return 1; -} - - -void MeshTopology :: GetFaceVertices (int fnr, Array & vertices) const -{ - vertices.SetSize(4); - int i; - for (i = 1; i <= 4; i++) - vertices.Elem(i) = face2vert.Get(fnr)[i-1]; - if (vertices.Elem(4) == 0) - vertices.SetSize(3); -} - -void MeshTopology :: GetFaceVertices (int fnr, int * vertices) const -{ - for (int i = 0; i <= 3; i++) - vertices[i] = face2vert.Get(fnr)[i]; -} - - -void MeshTopology :: GetEdgeVertices (int ednr, int & v1, int & v2) const -{ - v1 = edge2vert.Get(ednr)[0]; - v2 = edge2vert.Get(ednr)[1]; -} - - -void MeshTopology :: GetFaceEdges (int fnr, Array & fedges, bool withorientation) const -{ - ArrayMem pi(4); - ArrayMem eledges; - - fedges.SetSize (0); - GetFaceVertices(fnr, pi); - - // Sort Edges according to global vertex numbers - // e1 = fmax, f2 - // e2 = fmax, f1 - // e3 = op e1(f2,f3) - // e4 = op e2(f1,f3) - - /* ArrayMem fp; - fp[0] = pi[0]; - for(int k=1;kfp[0]) swap(fp[k],fp[0]); - - fp[1] = fp[0]+ */ - - - // GetVertexElements (pi[0], els); - FlatArray els= GetVertexElements (pi[0]); - - // find one element having all vertices of the face - for (int i = 0; i < els.Size(); i++) - { - const Element & el = mesh.VolumeElement(els[i]); - int nref_faces = GetNFaces (el.GetType()); - const ELEMENT_FACE * ref_faces = GetFaces (el.GetType()); - int nfa_ref_edges = GetNEdges (GetFaceType(fnr)); - - int cntv = 0,fa=-1; - for(int m=0;m0;j++) - for(int k=0;k 0 ) ? 1 : -1; } - if (cntv == pi.Size()) - { - fa=m; - break; - } - } + } + else + { + for (i = 0; i < 4; i++) + { + if (!surfedges.Get(elnr)[i]) return i; + eledges[i] = abs (surfedges.Get(elnr)[i]); + } + } + return 4; + } + else + { + eledges[0] = abs (segedges.Get(elnr)); + if (orient) + orient[0] = segedges.Get(elnr) > 0 ? 1 : -1; + } + return 1; + } + + + void MeshTopology :: GetFaceVertices (int fnr, Array & vertices) const + { + vertices.SetSize(4); + int i; + for (i = 1; i <= 4; i++) + vertices.Elem(i) = face2vert.Get(fnr)[i-1]; + if (vertices.Elem(4) == 0) + vertices.SetSize(3); + } + + void MeshTopology :: GetFaceVertices (int fnr, int * vertices) const + { + for (int i = 0; i <= 3; i++) + vertices[i] = face2vert.Get(fnr)[i]; + } + + + void MeshTopology :: GetEdgeVertices (int ednr, int & v1, int & v2) const + { + v1 = edge2vert.Get(ednr)[0]; + v2 = edge2vert.Get(ednr)[1]; + } + + + void MeshTopology :: GetFaceEdges (int fnr, Array & fedges, bool withorientation) const + { + ArrayMem pi(4); + ArrayMem eledges; + + fedges.SetSize (0); + GetFaceVertices(fnr, pi); + + // Sort Edges according to global vertex numbers + // e1 = fmax, f2 + // e2 = fmax, f1 + // e3 = op e1(f2,f3) + // e4 = op e2(f1,f3) + + /* ArrayMem fp; + fp[0] = pi[0]; + for(int k=1;kfp[0]) swap(fp[k],fp[0]); + + fp[1] = fp[0]+ */ + + + // GetVertexElements (pi[0], els); + FlatArray els= GetVertexElements (pi[0]); + + // find one element having all vertices of the face + for (int i = 0; i < els.Size(); i++) + { + const Element & el = mesh.VolumeElement(els[i]); + int nref_faces = GetNFaces (el.GetType()); + const ELEMENT_FACE * ref_faces = GetFaces1 (el.GetType()); + int nfa_ref_edges = GetNEdges (GetFaceType(fnr)); + + int cntv = 0,fa=-1; + for(int m=0;m0;j++) + for(int k=0;k=0) - { - const ELEMENT_EDGE * fa_ref_edges = GetEdges(GetFaceType(fnr)); - fedges.SetSize(nfa_ref_edges); - GetElementEdges (els[i], eledges); + if(fa>=0) + { + const ELEMENT_EDGE * fa_ref_edges = GetEdges1 (GetFaceType(fnr)); + fedges.SetSize(nfa_ref_edges); + GetElementEdges (els[i], eledges); - for (int j = 0; j < eledges.Size(); j++) - { - int vi1, vi2; - GetEdgeVertices (eledges[j], vi1, vi2); + for (int j = 0; j < eledges.Size(); j++) + { + int vi1, vi2; + GetEdgeVertices (eledges[j], vi1, vi2); - bool has1 = 0; - bool has2 = 0; - for (int k = 0; k < pi.Size(); k++) - { - if (vi1 == pi[k]) has1 = 1; - if (vi2 == pi[k]) has2 = 1; + bool has1 = 0; + bool has2 = 0; + for (int k = 0; k < pi.Size(); k++) + { + if (vi1 == pi[k]) has1 = 1; + if (vi2 == pi[k]) has2 = 1; - } + } - if (has1 && has2) // eledges[j] is on face - { - // fedges.Append (eledges[j]); - for(int k=0;k & elements) const -{ - if (vert2element) - { - int i; - int ne = vert2element->EntrySize(vnr); - elements.SetSize(ne); - for (i = 1; i <= ne; i++) - elements.Elem(i) = vert2element->Get(vnr, i); - } -} - - -FlatArray MeshTopology :: GetVertexElements (int vnr) const -{ - if (vert2element) - return (*vert2element)[vnr]; - return FlatArray (0,0); -} - -FlatArray MeshTopology :: GetVertexSurfaceElements (int vnr) const -{ - if (vert2surfelement) - return (*vert2surfelement)[vnr]; - return FlatArray (0,0); -} - - -void MeshTopology :: GetVertexSurfaceElements( int vnr, - Array& elements ) const -{ - if (vert2surfelement) - { - int i; - int ne = vert2surfelement->EntrySize(vnr); - elements.SetSize(ne); - for (i = 1; i <= ne; i++) - elements.Elem(i) = vert2surfelement->Get(vnr, i); - } -} - - -int MeshTopology :: GetVerticesEdge ( int v1, int v2 ) const -{ - Array elements_v1, elementedges; - GetVertexElements ( v1, elements_v1); - int edv1, edv2; - - for ( int i = 0; i < elements_v1.Size(); i++ ) - { - GetElementEdges( elements_v1[i], elementedges ); - for ( int ed = 0; ed < elementedges.Size(); ed ++) - { - GetEdgeVertices( elementedges[ed], edv1, edv2 ); - if ( ( edv1 == v1 && edv2 == v2 ) || ( edv1 == v2 && edv2 == v1 ) ) - return elementedges[ed]; - } - } - - return -1; -} - - - -void MeshTopology :: GetSegmentVolumeElements ( int segnr, Array & volels ) const -{ - int v1, v2; - GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 ); - Array volels1, volels2; - GetVertexElements ( v1, volels1 ); - GetVertexElements ( v2, volels2 ); - volels.SetSize(0); - - for ( int eli1=1; eli1 <= volels1.Size(); eli1++) - if ( volels2.Contains( volels1.Elem(eli1) ) ) - volels.Append ( volels1.Elem(eli1) ); - -} + return; + } + } + } + + + ELEMENT_TYPE MeshTopology :: GetFaceType (int fnr) const + { + if (face2vert.Get(fnr)[3] == 0) return TRIG; else return QUAD; + } + + + void MeshTopology :: GetVertexElements (int vnr, Array & elements) const + { + if (vert2element) + { + int i; + int ne = vert2element->EntrySize(vnr); + elements.SetSize(ne); + for (i = 1; i <= ne; i++) + elements.Elem(i) = vert2element->Get(vnr, i); + } + } + + + FlatArray MeshTopology :: GetVertexElements (int vnr) const + { + if (vert2element) + return (*vert2element)[vnr]; + return FlatArray (0,0); + } + + FlatArray MeshTopology :: GetVertexSurfaceElements (int vnr) const + { + if (vert2surfelement) + return (*vert2surfelement)[vnr]; + return FlatArray (0,0); + } + + + void MeshTopology :: GetVertexSurfaceElements( int vnr, + Array& elements ) const + { + if (vert2surfelement) + { + int i; + int ne = vert2surfelement->EntrySize(vnr); + elements.SetSize(ne); + for (i = 1; i <= ne; i++) + elements.Elem(i) = vert2surfelement->Get(vnr, i); + } + } + + + int MeshTopology :: GetVerticesEdge ( int v1, int v2 ) const + { + Array elements_v1, elementedges; + GetVertexElements ( v1, elements_v1); + int edv1, edv2; + + for ( int i = 0; i < elements_v1.Size(); i++ ) + { + GetElementEdges( elements_v1[i], elementedges ); + for ( int ed = 0; ed < elementedges.Size(); ed ++) + { + GetEdgeVertices( elementedges[ed], edv1, edv2 ); + if ( ( edv1 == v1 && edv2 == v2 ) || ( edv1 == v2 && edv2 == v1 ) ) + return elementedges[ed]; + } + } + + return -1; + } + + + + void MeshTopology :: GetSegmentVolumeElements ( int segnr, Array & volels ) const + { + int v1, v2; + GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 ); + Array volels1, volels2; + GetVertexElements ( v1, volels1 ); + GetVertexElements ( v2, volels2 ); + volels.SetSize(0); + + for ( int eli1=1; eli1 <= volels1.Size(); eli1++) + if ( volels2.Contains( volels1.Elem(eli1) ) ) + volels.Append ( volels1.Elem(eli1) ); + + } } diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 6a8d321c..d46342a2 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -62,8 +62,10 @@ public: static inline int GetNFaces (ELEMENT_TYPE et); static const Point3d * GetVertices (ELEMENT_TYPE et); - inline static const ELEMENT_EDGE * GetEdges (ELEMENT_TYPE et); - inline static const ELEMENT_FACE * GetFaces (ELEMENT_TYPE et); + inline static const ELEMENT_EDGE * GetEdges1 (ELEMENT_TYPE et); + inline static const ELEMENT_EDGE * GetEdges0 (ELEMENT_TYPE et); + inline static const ELEMENT_FACE * GetFaces1 (ELEMENT_TYPE et); + inline static const ELEMENT_FACE * GetFaces0 (ELEMENT_TYPE et); int GetSegmentEdge (int segnr) const { return abs(segedges[segnr-1]); } @@ -295,7 +297,7 @@ int MeshTopology :: GetNFaces (ELEMENT_TYPE et) -const ELEMENT_EDGE * MeshTopology :: GetEdges (ELEMENT_TYPE et) +const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et) { static int segm_edges[1][2] = { { 1, 2 }}; @@ -392,7 +394,113 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges (ELEMENT_TYPE et) } -const ELEMENT_FACE * MeshTopology :: GetFaces (ELEMENT_TYPE et) + +const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) +{ + static int segm_edges[1][2] = + { { 0, 1 }}; + + static int trig_edges[3][2] = + { { 2, 0 }, + { 1, 2 }, + { 0, 1 }}; + + static int quad_edges[4][2] = + { { 0, 1 }, + { 2, 3 }, + { 3, 0 }, + { 1, 2 }}; + + + static int tet_edges[6][2] = + { { 3, 0 }, + { 3, 1 }, + { 3, 2 }, + { 0, 1 }, + { 0, 2 }, + { 1, 2 }}; + + static int prism_edges[9][2] = + { { 2, 0 }, + { 0, 1 }, + { 2, 1 }, + { 5, 3 }, + { 3, 4 }, + { 5, 4 }, + { 2, 5 }, + { 0, 3 }, + { 1, 4 }}; + + static int pyramid_edges[8][2] = + { { 0, 1 }, + { 1, 2 }, + { 0, 3 }, + { 3, 2 }, + { 0, 4 }, + { 1, 4 }, + { 2, 4 }, + { 3, 4 }}; + + static int hex_edges[12][2] = + { + { 0, 1 }, + { 2, 3 }, + { 3, 0 }, + { 1, 2 }, + { 4, 5 }, + { 6, 7 }, + { 7, 4 }, + { 5, 6 }, + { 0, 4 }, + { 1, 5 }, + { 2, 6 }, + { 3, 7 }, + }; + + switch (et) + { + case SEGMENT: + case SEGMENT3: + return segm_edges; + + case TRIG: + case TRIG6: + return trig_edges; + + case QUAD: + case QUAD6: + case QUAD8: + return quad_edges; + + case TET: + case TET10: + return tet_edges; + + case PYRAMID: + return pyramid_edges; + + case PRISM: + case PRISM12: + return prism_edges; + + case HEX: + return hex_edges; + default: + cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; + } + return 0; +} + + + + + + + + + + +const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) { static const int trig_faces[1][4] = { { 1, 2, 3, 0 } }; @@ -474,6 +582,89 @@ const ELEMENT_FACE * MeshTopology :: GetFaces (ELEMENT_TYPE et) +const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) +{ + static const int trig_faces[1][4] = + { { 0, 1, 2, -1 } }; + static const int quad_faces[1][4] = + { { 0, 1, 2, 3 } }; + + static const int tet_faces[4][4] = + { { 3, 1, 2, -1 }, + { 3, 2, 0, -1 }, + { 3, 0, 1, -1 }, + { 0, 2, 1, -1 } }; + + static const int prism_faces[5][4] = + { + { 0, 2, 1, -1 }, + { 3, 4, 5, -1 }, + { 2, 0, 3, 5 }, + { 0, 1, 4, 3 }, + { 1, 2, 5, 4 } + }; + + static const int pyramid_faces[5][4] = + { + { 0, 1, 4, -1 }, + { 1, 2, 4, -1 }, + { 2, 3, 4, -1 }, + { 3, 0, 4, -1 }, + { 0, 3, 2, 1 } + }; + + static const int hex_faces[6][4] = + { + { 0, 3, 2, 1 }, + { 4, 5, 6, 7 }, + { 0, 1, 5, 4 }, + { 1, 2, 6, 5 }, + { 2, 3, 7, 6 }, + { 3, 0, 4, 7 } + }; + + + + switch (et) + { + case TRIG: + case TRIG6: + return trig_faces; + + case QUAD: + case QUAD6: + case QUAD8: + return quad_faces; + + + case TET: + case TET10: + return tet_faces; + + case PRISM: + case PRISM12: + return prism_faces; + + case PYRAMID: + return pyramid_faces; + + case SEGMENT: + case SEGMENT3: + + case HEX: + return hex_faces; + + default: + cerr << "Ng_ME_GetVertices, illegal element type " << et << endl; + } + return 0; +} + + + + + + diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index d7b13338..39bf7458 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -1880,7 +1880,7 @@ namespace netgen if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) { - const ELEMENT_FACE * faces = MeshTopology :: GetFaces (TET); + const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (TET); const Point3d * vertices = MeshTopology :: GetVertices (TET); /* @@ -2084,7 +2084,7 @@ namespace netgen CurvedElements & curv = mesh->GetCurvedElements(); if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) { - const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PRISM); + const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (PRISM); const Point3d * vertices = MeshTopology :: GetVertices (PRISM); Point<3> grid[11][11]; @@ -2469,7 +2469,7 @@ namespace netgen glEnd (); */ - const ELEMENT_FACE * faces = MeshTopology :: GetFaces (HEX); + const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (HEX); const Point3d * vertices = MeshTopology :: GetVertices (HEX); Point<3> grid[11][11]; @@ -2618,7 +2618,7 @@ namespace netgen if (curv.IsHighOrder()) // && curv.IsElementCurved(ei)) { - const ELEMENT_FACE * faces = MeshTopology :: GetFaces (PYRAMID); + const ELEMENT_FACE * faces = MeshTopology :: GetFaces1 (PYRAMID); const Point3d * vertices = MeshTopology :: GetVertices (PYRAMID); Point<3> grid[11][11];