export Trafo to py, IdentifyPeriodic with Trafo

This commit is contained in:
Joachim Schöberl 2018-03-09 15:29:50 +01:00
parent ee7ac2e0a0
commit aa8dbac6be
5 changed files with 37 additions and 14 deletions

View File

@ -120,9 +120,11 @@ PeriodicIdentification ::
PeriodicIdentification (int anr,
const CSGeometry & ageom,
const Surface * as1,
const Surface * as2)
: Identification(anr, ageom)
const Surface * as2,
Transformation<3> atrafo)
: Identification(anr, ageom), trafo(atrafo)
{
inv_trafo = trafo.CalcInverse();
s1 = as1;
s2 = as2;
}
@ -256,25 +258,26 @@ GetIdentifiedPoint (class Mesh & mesh, int pi)
const Surface *snew;
const Point<3> & p = mesh.Point (pi);
Point<3> hp = p;
if (s1->PointOnSurface (p))
{
snew = s2;
hp = trafo(hp);
}
else
{
if (s2->PointOnSurface (p))
{
snew = s1;
hp = inv_trafo(hp);
}
else
{
cerr << "GetIdenfifiedPoint: Not possible" << endl;
exit (1);
throw NgException("GetIdenfifiedPoint: Not possible");
}
}
// project to other surface
Point<3> hp = p;
snew->Project (hp);
int newpi = 0;
@ -306,15 +309,15 @@ GetIdentifiedPoint (class Mesh & mesh, int pi)
void PeriodicIdentification :: IdentifyPoints (class Mesh & mesh)
{
int i, j;
for (i = 1; i <= mesh.GetNP(); i++)
for (int i = 1; i <= mesh.GetNP(); i++)
{
Point<3> p = mesh.Point(i);
if (s1->PointOnSurface (p))
{
Point<3> pp = p;
pp = trafo(pp);
s2->Project (pp);
for (j = 1; j <= mesh.GetNP(); j++)
for (int j = 1; j <= mesh.GetNP(); j++)
if (Dist2(mesh.Point(j), pp) < 1e-6)
{
mesh.GetIdentifications().Add (i, j, nr);

View File

@ -78,11 +78,14 @@ namespace netgen
{
const Surface * s1;
const Surface * s2;
Transformation<3> trafo; // from s1 to s2
Transformation<3> inv_trafo; // from s2 to s1
public:
PeriodicIdentification (int anr,
const CSGeometry & ageom,
const Surface * as1,
const Surface * as2);
const Surface * as2,
Transformation<3> atrafo = Vec<3>(0,0,0));
virtual ~PeriodicIdentification ();
virtual void Print (ostream & ost) const;
virtual void GetData (ostream & ost) const;

View File

@ -555,7 +555,8 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
)
.def("PeriodicSurfaces", FunctionPointer
([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2,
Transformation<3> trafo)
{
Array<int> si1, si2;
s1->GetSolid()->GetSurfaceIndices (si1);
@ -564,9 +565,11 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails!
self.AddIdentification
(new PeriodicIdentification
(self.GetNIdentifications()+1, self,
self.GetSurface (si1[0]), self.GetSurface (si2[0])));
self.GetSurface (si1[0]), self.GetSurface (si2[0]),
trafo));
}),
py::arg("solid1"), py::arg("solid2")
py::arg("solid1"), py::arg("solid2"),
py::arg("trafo")=Transformation<3>(Vec<3>(0,0,0))
)
.def("AddPoint", [] (CSGeometry & self, Point<3> p, int index) -> CSGeometry&

View File

@ -129,7 +129,14 @@ public:
}
///
void CalcInverse (Transformation & inv) const;
Transformation CalcInverse () const
{
Transformation inv;
// inv.m = Inv(m);
::netgen::CalcInverse (m, inv.m);
inv.v = inv.m * (-v);
return inv;
}
/// this = ta x tb
void Combine (const Transformation & ta, const Transformation & tb)

View File

@ -86,11 +86,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
m.def ("Pnt", FunctionPointer
([](double x, double y) { return Point<2>(x,y); }));
/*
// duplicated functions ????
m.def ("Pnt", FunctionPointer
([](double x, double y, double z) { return Point<3>(x,y,z); }));
m.def ("Pnt", FunctionPointer
([](double x, double y) { return Point<2>(x,y); }));
*/
py::class_<Vec<2>> (m, "Vec2d")
.def(py::init<double,double>())
@ -117,6 +119,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
m.def ("Vec", FunctionPointer
([] (double x, double y) { return Vec<2>(x,y); }));
py::class_<Transformation<3>> (m, "Trafo")
.def(py::init<Vec<3>>())
.def("__call__", [] (Transformation<3> trafo, Point<3> p) { return trafo(p); })
// .def ("__str__", &ToString<Transformation<3>>)
;
m.def ("SetTransformation", FunctionPointer
([](int dir, double angle)