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

mesh identify periodic for non tet meshes

Browse Source
This commit is contained in:
Christopher Lackner 2020-03-19 18:12:55 +01:00
parent bff0e67576
commit 1f78f900dd
3 changed files with 19 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
libsrc/meshing/meshclass.cpp
View File

@ -6115,15 +6115,19 @@ namespace netgen
int Mesh::IdentifyPeriodicBoundaries(const string &s1, int Mesh::IdentifyPeriodicBoundaries(const string &s1,
const string &s2, const string &s2,
const Transformation<3> &mapping) const Transformation<3> &mapping,
double pointTolerance)
{ {
auto nr = ident->GetMaxNr() + 1; auto nr = ident->GetMaxNr() + 1;
ident->SetType(nr, Identifications::PERIODIC); ident->SetType(nr, Identifications::PERIODIC);
double lami[4]; double lami[4];
set<int> identified_points; set<int> identified_points;
Point3d pmin, pmax; if(pointTolerance < 0.)
GetBox(pmin, pmax); {
auto eps = 1e-8 * (pmax-pmin).Length(); Point3d pmin, pmax;
GetBox(pmin, pmax);
pointTolerance = 1e-8 * (pmax-pmin).Length();
}
for(const auto& se : surfelements) for(const auto& se : surfelements)
{ {
if(GetBCName(se.index-1) != s1) if(GetBCName(se.index-1) != s1)
@ -6138,14 +6142,18 @@ namespace netgen
auto other_nr = GetElementOfPoint(mapped_pt, lami, true); auto other_nr = GetElementOfPoint(mapped_pt, lami, true);
int index = -1; int index = -1;
auto other_el = VolumeElement(other_nr); auto other_el = VolumeElement(other_nr);
for(auto i : Range(4)) for(auto i : Range(other_el.PNums().Size()))
if((mapped_pt - (*this)[other_el.PNums()[i]]).Length() < eps) if((mapped_pt - (*this)[other_el.PNums()[i]]).Length() < pointTolerance)
{ {
index = i; index = i;
break; break;
} }
if(index == -1) if(index == -1)
throw Exception("Did not find mapped point, are you sure your mesh is periodic?"); {
cout << "point coordinates = " << pt << endl;
cout << "mapped coordinates = " << mapped_pt << endl;
throw Exception("Did not find mapped point with nr " + ToString(pi) + ", are you sure your mesh is periodic?");
}
auto other_pi = other_el.PNums()[index]; auto other_pi = other_el.PNums()[index];
identified_points.insert(pi); identified_points.insert(pi);
ident->Add(pi, other_pi, nr); ident->Add(pi, other_pi, nr);

3
libsrc/meshing/meshclass.hpp
View File

@ -712,7 +712,8 @@ namespace netgen
int IdentifyPeriodicBoundaries(const string& s1, int IdentifyPeriodicBoundaries(const string& s1,
const string& s2, const string& s2,
const Transformation<3>& mapping); const Transformation<3>& mapping,
double pointTolerance);
// #ifdef NONE // #ifdef NONE
// /* // /*

3
libsrc/meshing/python_mesh.cpp
View File

@ -855,7 +855,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
py::arg("pid2"), py::arg("pid2"),
py::arg("identnr"), py::arg("identnr"),
py::arg("type")) py::arg("type"))
.def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries) .def("IdentifyPeriodicBoundaries", &Mesh::IdentifyPeriodicBoundaries,
py::arg("face1"), py::arg("face2"), py::arg("mapping"), py::arg("point_tolerance") = -1.)
.def("GetNrIdentifications", [](Mesh& self) .def("GetNrIdentifications", [](Mesh& self)
{ {
return self.GetIdentifications().GetMaxNr(); return self.GetIdentifications().GetMaxNr();