From 4b26f399741acf2f4dcd77c9bbf10249b13f4525 Mon Sep 17 00:00:00 2001 From: Michael Neunteufel Date: Fri, 13 Sep 2019 09:09:56 +0200 Subject: [PATCH 1/2] export hprefleft/hprefright to python, sort points for segments and check if point is part of correct material --- libsrc/geom2d/python_geom2d.cpp | 8 ++++---- libsrc/meshing/classifyhpel.hpp | 29 ++++++++++++++++++----------- libsrc/meshing/hprefinement.cpp | 33 ++++++++++++++++++++------------- 3 files changed, 42 insertions(+), 28 deletions(-) diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index e2089d3e..54f28ef1 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -63,7 +63,7 @@ DLL_HEADER void ExportGeom2d(py::module &m) py::arg("x"), py::arg("y"), py::arg("maxh") = 1e99, py::arg("hpref")=0, py::arg("name")="") .def("Append", FunctionPointer([](SplineGeometry2d &self, py::list segment, int leftdomain, int rightdomain, optional> bc, optional copy, double maxh, - double hpref) + double hpref, double hprefleft, double hprefright) { auto segtype = py::cast(segment[0]); @@ -86,8 +86,8 @@ DLL_HEADER void ExportGeom2d(py::module &m) seg->leftdom = leftdomain; seg->rightdom = rightdomain; seg->hmax = maxh; - seg->hpref_left = hpref; - seg->hpref_right = hpref; + seg->hpref_left = max(hpref, hprefleft); + seg->hpref_right = max(hpref,hprefright); seg->reffak = 1; seg->copyfrom = -1; if (copy.has_value()) @@ -110,7 +110,7 @@ DLL_HEADER void ExportGeom2d(py::module &m) return self.GetNSplines()-1; }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = py::int_(0), py::arg("bc")=nullopt, py::arg("copy")=nullopt, py::arg("maxh")=1e99, - py::arg("hpref")=0) + py::arg("hpref")=0,py::arg("hprefleft")=0,py::arg("hprefright")=0) .def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, py::list point_indices, int leftdomain, int rightdomain) diff --git a/libsrc/meshing/classifyhpel.hpp b/libsrc/meshing/classifyhpel.hpp index 9eaea354..a44ad2b3 100644 --- a/libsrc/meshing/classifyhpel.hpp +++ b/libsrc/meshing/classifyhpel.hpp @@ -664,8 +664,14 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int dim, const FaceDescriptor & fd) { + cout << "IN ClassifyTrig!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; HPREF_ELEMENT_TYPE type = HP_NONE; - + cout << "el.index = " << el.index << endl; + cout << "edges = " << edges << endl; + cout << "edgepoint_dom = " << edgepoint_dom << endl; + cout << "face_edges = " << face_edges << endl; + cout << "surf_edges = " << surf_edges << endl; + cout << "facepoint = " << facepoint << endl; int pnums[3]; int p[3]; @@ -686,9 +692,9 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge { p[m] = (j+m)%3 +1; // local vertex number pnums[m] = el.PNum(p[m]); // global vertex number - // *testout << pnums[m] << " \t "; + cout << pnums[m] << " \t "; } - // *testout << endl ; + cout << endl ; if(dim == 3) { @@ -761,11 +767,12 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge int ep1=p[eledges[k][0]-1]; int ep2=p[eledges[k][1]-1]; - INDEX_2 i2(el.PNum(ep1),el.PNum(ep2)); + INDEX_2 i2 = INDEX_2::Sort(el.PNum(ep1),el.PNum(ep2)); + cout << "ep1 = " << ep1 << ", ep2 = " << ep2 << endl; if(edges.Used(i2)) { - + if(edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep1-1])) || edgepoint_dom.Used(INDEX_2(-1,pnums[ep1-1])) || edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep2-1])) || @@ -783,11 +790,10 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge for(int k=0;k<3;k++) - if(edgepoint.Test(pnums[k])) //edgepoint, but not member of sing_edge on trig -> cp + if(edgepoint.Test(pnums[k]) && (edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[k])) || edgepoint_dom.Used(INDEX_2(-1,pnums[k])))) //edgepoint, but not member of sing_edge on trig -> cp { INDEX_2 i2a=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+1)%3])); - INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3])); - + INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3])); if(!edges.Used(i2a) && !edges.Used(i2b)) point_sing[p[k]-1] = 3; } @@ -796,8 +802,8 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge if(cornerpoint.Test(el.PNum(p[k]))) point_sing[p[k]-1] = 3; - *testout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl; - + cout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl; + cout << "edge_sing = " << edge_sing[0] << edge_sing[1] << edge_sing[2] << endl; if(edge_sing[0] + edge_sing[1] + edge_sing[2] == 0) { int ps = point_sing[0] + point_sing[1] + point_sing[2]; @@ -860,7 +866,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge if(type!=HP_NONE) break; } - *testout << "type = " << type << endl; + cout << "type = " << type << endl; for(int k=0;k<3;k++) el[k] = pnums[k]; /*if(type != HP_NONE) @@ -871,6 +877,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge cout << " type " << type << endl; } */ + cout << "End ClassifyTrig!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; return(type); } #ifdef HPREF_OLD diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 8d046c01..5a31b223 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -615,6 +615,7 @@ namespace netgen void DoRefinement (Mesh & mesh, NgArray & elements, Refinement * ref, double fac1) { + cout << "IN DOREFINEMENT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; elements.SetAllocSize (5 * elements.Size()); INDEX_2_HASHTABLE newpts(elements.Size()+1); INDEX_3_HASHTABLE newfacepts(elements.Size()+1); @@ -1303,6 +1304,7 @@ namespace netgen void HPRefinement (Mesh & mesh, Refinement * ref, int levels, double fac1, bool setorders, bool reflevels) { + cout << "IN HPRefinement!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; PrintMessage (1, "HP Refinement called, levels = ", levels); @@ -1560,7 +1562,8 @@ namespace netgen bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HASHTABLE & edgepoint_dom, NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE & faces, INDEX_2_HASHTABLE & face_edges, INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int & levels, int & act_ref) -{ +{ + cout << "IN CheckSingularities!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; bool sing = 0; if (mesh.GetDimension() == 3) { @@ -1651,8 +1654,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS } edgepoint.Or (cornerpoint); - (*testout) << "cornerpoint = " << endl << cornerpoint << endl; - (*testout) << "edgepoint = " << endl << edgepoint << endl; + cout << "cornerpoint = " << endl << cornerpoint << endl; + cout << "edgepoint = " << endl << edgepoint << endl; facepoint = 0; for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) @@ -1722,14 +1725,14 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (seg.singedge_left * levels >= act_ref) { - INDEX_2 i2 (mesh.LineSegment(i)[0], + INDEX_2 i2 = INDEX_2::Sort(mesh.LineSegment(i)[0], mesh.LineSegment(i)[1]); edges.Set(i2,1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); - *testout << " singleft " << endl; - *testout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; - *testout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; + cout << " singleft " << endl; + cout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; + cout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domin, i2.I1()), 1); edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domin, i2.I2()), 1); sing = 1; @@ -1738,15 +1741,15 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (seg.singedge_right * levels >= act_ref) { - INDEX_2 i2 (mesh.LineSegment(i)[1], + INDEX_2 i2 = INDEX_2::Sort(mesh.LineSegment(i)[1], mesh.LineSegment(i)[0]); edges.Set (i2, 1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); - *testout << " singright " << endl; - *testout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; - *testout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; + cout << " singright " << endl; + cout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; + cout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domout, i2.I1()), 1); edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domout, i2.I2()), 1); @@ -1794,8 +1797,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS edgepoint.Or (cornerpoint); - (*testout) << "2d sing edges: " << endl << edges << endl; - (*testout) << "2d cornerpoints: " << endl << cornerpoint << endl + cout << "2d sing edges: " << endl << edges << endl; + cout << "2d cornerpoints: " << endl << cornerpoint << endl << "2d edgepoints: " << endl << edgepoint << endl; facepoint = 0; @@ -1804,6 +1807,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (!sing) cout << "PrepareElements no more to do for actual refinement " << act_ref << endl; + cout << "End CheckSingularities!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; + return(sing); } @@ -1811,6 +1816,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS bool ClassifyHPElements (Mesh & mesh, NgArray & elements, int & act_ref, int & levels) { + cout << "IN ClassifyHPElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; INDEX_2_HASHTABLE edges(mesh.GetNSeg()+1); NgBitArray edgepoint(mesh.GetNP()); INDEX_2_HASHTABLE edgepoint_dom(mesh.GetNSeg()+1); @@ -1956,6 +1962,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (misses[i]) cout << " in update classification missing case " << i << " occurred " << misses[i] << " times" << endl; + cout << "end ClassifyHPElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; return(sing); } } From 9f79451fe1b8273d39c5689b7f85fb4ac1ad8a3a Mon Sep 17 00:00:00 2001 From: Michael Neunteufel Date: Fri, 13 Sep 2019 09:42:05 +0200 Subject: [PATCH 2/2] clean up --- libsrc/meshing/classifyhpel.hpp | 28 ++++++++++------------------ libsrc/meshing/hprefinement.cpp | 27 ++++++++++----------------- 2 files changed, 20 insertions(+), 35 deletions(-) diff --git a/libsrc/meshing/classifyhpel.hpp b/libsrc/meshing/classifyhpel.hpp index a44ad2b3..259a7cd0 100644 --- a/libsrc/meshing/classifyhpel.hpp +++ b/libsrc/meshing/classifyhpel.hpp @@ -664,14 +664,8 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int dim, const FaceDescriptor & fd) { - cout << "IN ClassifyTrig!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; - HPREF_ELEMENT_TYPE type = HP_NONE; - cout << "el.index = " << el.index << endl; - cout << "edges = " << edges << endl; - cout << "edgepoint_dom = " << edgepoint_dom << endl; - cout << "face_edges = " << face_edges << endl; - cout << "surf_edges = " << surf_edges << endl; - cout << "facepoint = " << facepoint << endl; + HPREF_ELEMENT_TYPE type = HP_NONE; + int pnums[3]; int p[3]; @@ -692,9 +686,9 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge { p[m] = (j+m)%3 +1; // local vertex number pnums[m] = el.PNum(p[m]); // global vertex number - cout << pnums[m] << " \t "; + // *testout << pnums[m] << " \t "; } - cout << endl ; + // *testout << endl ; if(dim == 3) { @@ -768,11 +762,9 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge int ep2=p[eledges[k][1]-1]; INDEX_2 i2 = INDEX_2::Sort(el.PNum(ep1),el.PNum(ep2)); - cout << "ep1 = " << ep1 << ", ep2 = " << ep2 << endl; if(edges.Used(i2)) { - if(edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep1-1])) || edgepoint_dom.Used(INDEX_2(-1,pnums[ep1-1])) || edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[ep2-1])) || @@ -793,17 +785,18 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge if(edgepoint.Test(pnums[k]) && (edgepoint_dom.Used(INDEX_2(fd.SurfNr(),pnums[k])) || edgepoint_dom.Used(INDEX_2(-1,pnums[k])))) //edgepoint, but not member of sing_edge on trig -> cp { INDEX_2 i2a=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+1)%3])); - INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3])); + INDEX_2 i2b=INDEX_2::Sort(el.PNum(p[k]), el.PNum(p[(k+2)%3])); + if(!edges.Used(i2a) && !edges.Used(i2b)) point_sing[p[k]-1] = 3; } for(int k=0;k<3;k++) if(cornerpoint.Test(el.PNum(p[k]))) - point_sing[p[k]-1] = 3; + point_sing[p[k]-1] = 3; - cout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl; - cout << "edge_sing = " << edge_sing[0] << edge_sing[1] << edge_sing[2] << endl; + *testout << "point_sing = " << point_sing[0] << point_sing[1] << point_sing[2] << endl; + if(edge_sing[0] + edge_sing[1] + edge_sing[2] == 0) { int ps = point_sing[0] + point_sing[1] + point_sing[2]; @@ -866,7 +859,7 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge if(type!=HP_NONE) break; } - cout << "type = " << type << endl; + *testout << "type = " << type << endl; for(int k=0;k<3;k++) el[k] = pnums[k]; /*if(type != HP_NONE) @@ -877,7 +870,6 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE & edge cout << " type " << type << endl; } */ - cout << "End ClassifyTrig!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; return(type); } #ifdef HPREF_OLD diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index 5a31b223..3aef7ebd 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -615,7 +615,6 @@ namespace netgen void DoRefinement (Mesh & mesh, NgArray & elements, Refinement * ref, double fac1) { - cout << "IN DOREFINEMENT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; elements.SetAllocSize (5 * elements.Size()); INDEX_2_HASHTABLE newpts(elements.Size()+1); INDEX_3_HASHTABLE newfacepts(elements.Size()+1); @@ -1304,7 +1303,6 @@ namespace netgen void HPRefinement (Mesh & mesh, Refinement * ref, int levels, double fac1, bool setorders, bool reflevels) { - cout << "IN HPRefinement!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; PrintMessage (1, "HP Refinement called, levels = ", levels); @@ -1563,7 +1561,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE & faces, INDEX_2_HASHTABLE & face_edges, INDEX_2_HASHTABLE & surf_edges, NgArray & facepoint, int & levels, int & act_ref) { - cout << "IN CheckSingularities!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; bool sing = 0; if (mesh.GetDimension() == 3) { @@ -1654,8 +1651,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS } edgepoint.Or (cornerpoint); - cout << "cornerpoint = " << endl << cornerpoint << endl; - cout << "edgepoint = " << endl << edgepoint << endl; + (*testout) << "cornerpoint = " << endl << cornerpoint << endl; + (*testout) << "edgepoint = " << endl << edgepoint << endl; facepoint = 0; for (SurfaceElementIndex sei = 0; sei < mesh.GetNSE(); sei++) @@ -1730,9 +1727,9 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS edges.Set(i2,1); edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); - cout << " singleft " << endl; - cout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; - cout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; + *testout << " singleft " << endl; + *testout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; + *testout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domin, i2.I1()), 1); edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domin, i2.I2()), 1); sing = 1; @@ -1747,9 +1744,9 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS edgepoint.Set(i2.I1()); edgepoint.Set(i2.I2()); - cout << " singright " << endl; - cout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; - cout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; + *testout << " singright " << endl; + *testout << " mesh.LineSegment(i).domout " << mesh.LineSegment(i).domout << endl; + *testout << " mesh.LineSegment(i).domin " << mesh.LineSegment(i).domin << endl; edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domout, i2.I1()), 1); edgepoint_dom.Set (INDEX_2(mesh.LineSegment(i).domout, i2.I2()), 1); @@ -1797,8 +1794,8 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS edgepoint.Or (cornerpoint); - cout << "2d sing edges: " << endl << edges << endl; - cout << "2d cornerpoints: " << endl << cornerpoint << endl + (*testout) << "2d sing edges: " << endl << edges << endl; + (*testout) << "2d cornerpoints: " << endl << cornerpoint << endl << "2d edgepoints: " << endl << edgepoint << endl; facepoint = 0; @@ -1807,8 +1804,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (!sing) cout << "PrepareElements no more to do for actual refinement " << act_ref << endl; - cout << "End CheckSingularities!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; - return(sing); } @@ -1816,7 +1811,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS bool ClassifyHPElements (Mesh & mesh, NgArray & elements, int & act_ref, int & levels) { - cout << "IN ClassifyHPElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; INDEX_2_HASHTABLE edges(mesh.GetNSeg()+1); NgBitArray edgepoint(mesh.GetNP()); INDEX_2_HASHTABLE edgepoint_dom(mesh.GetNSeg()+1); @@ -1962,7 +1956,6 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE & edges, INDEX_2_HAS if (misses[i]) cout << " in update classification missing case " << i << " occurred " << misses[i] << " times" << endl; - cout << "end ClassifyHPElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl; return(sing); } }