ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2024-11-11 16:49:16 +05:00
Code Issues Packages Projects Releases Wiki Activity

export hprefleft/hprefright to python, sort points for segments and check if point is part of correct material

Browse Source
This commit is contained in:
Michael Neunteufel 2019-09-13 09:09:56 +02:00
parent c359832767
commit 4b26f39974
3 changed files with 42 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
libsrc/geom2d/python_geom2d.cpp
View File

@ -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<variant<int, string>> bc, optional<int> copy, double maxh,
double hpref)
double hpref, double hprefleft, double hprefright)
{
auto segtype = py::cast<std::string>(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)

29
libsrc/meshing/classifyhpel.hpp
View File

@ -664,8 +664,14 @@ HPREF_ELEMENT_TYPE ClassifyTrig(HPRefElement & el, INDEX_2_HASHTABLE<int> & edge
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & 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<int> & 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<int> & 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<int> & 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<int> & 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<int> & 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<int> & edge
cout << " type " << type << endl;
}
*/
cout << "End ClassifyTrig!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
return(type);
}
#ifdef HPREF_OLD

33
libsrc/meshing/hprefinement.cpp
View File

@ -615,6 +615,7 @@ namespace netgen
void DoRefinement (Mesh & mesh, NgArray<HPRefElement> & elements,
Refinement * ref, double fac1)
{
cout << "IN DOREFINEMENT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
elements.SetAllocSize (5 * elements.Size());
INDEX_2_HASHTABLE<int> newpts(elements.Size()+1);
INDEX_3_HASHTABLE<int> 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<int> & edges, INDEX_2_HASHTABLE<int> & edgepoint_dom,
NgBitArray & cornerpoint, NgBitArray & edgepoint, INDEX_3_HASHTABLE<int> & faces, INDEX_2_HASHTABLE<int> & face_edges,
INDEX_2_HASHTABLE<int> & surf_edges, NgArray<int, PointIndex::BASE> & 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<int> & 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<int> & 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<int> & 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<int> & 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<int> & 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<int> & edges, INDEX_2_HAS
bool ClassifyHPElements (Mesh & mesh, NgArray<HPRefElement> & elements, int & act_ref, int & levels)
{
cout << "IN ClassifyHPElements!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
INDEX_2_HASHTABLE<int> edges(mesh.GetNSeg()+1);
NgBitArray edgepoint(mesh.GetNP());
INDEX_2_HASHTABLE<int> edgepoint_dom(mesh.GetNSeg()+1);
@ -1956,6 +1962,7 @@ bool CheckSingularities(Mesh & mesh, INDEX_2_HASHTABLE<int> & 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);
}
}