From a4fe0c1c4168109bd34b6adbfd9144733a650a84 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 4 Oct 2016 19:30:57 +0200 Subject: [PATCH 1/9] first push --- libsrc/csg/CMakeLists.txt | 2 +- libsrc/csg/csg.hpp | 2 +- libsrc/csg/genmesh.cpp | 14 ++++ libsrc/csg/python_csg.cpp | 49 +++++++++++- libsrc/csg/splinesurface.cpp | 115 ++++++++++++++++++++++++++++ libsrc/csg/splinesurface.cpp~ | 92 ++++++++++++++++++++++ libsrc/csg/splinesurface.hpp | 47 ++++++++++++ libsrc/csg/splinesurface.hpp~ | 46 +++++++++++ libsrc/include/nginterface.h | 3 + libsrc/interface/nginterface.cpp | 5 ++ libsrc/interface/nginterface_v2.cpp | 16 ++++ libsrc/meshing/meshclass.cpp | 92 +++++++++++++++++++++- libsrc/meshing/meshclass.hpp | 8 ++ libsrc/meshing/meshtype.hpp | 4 +- 14 files changed, 488 insertions(+), 7 deletions(-) create mode 100644 libsrc/csg/splinesurface.cpp create mode 100644 libsrc/csg/splinesurface.cpp~ create mode 100644 libsrc/csg/splinesurface.hpp create mode 100644 libsrc/csg/splinesurface.hpp~ diff --git a/libsrc/csg/CMakeLists.txt b/libsrc/csg/CMakeLists.txt index 110c8bad..f1ff6c9c 100644 --- a/libsrc/csg/CMakeLists.txt +++ b/libsrc/csg/CMakeLists.txt @@ -5,7 +5,7 @@ add_library(csg ${NG_LIB_TYPE} explicitcurve2d.cpp extrusion.cpp gencyl.cpp genmesh.cpp identify.cpp manifold.cpp meshsurf.cpp polyhedra.cpp revolution.cpp singularref.cpp solid.cpp specpoin.cpp spline3d.cpp surface.cpp triapprox.cpp zrefine.cpp - python_csg.cpp + python_csg.cpp splinesurface.cpp ) if(APPLE) set_target_properties( csg PROPERTIES SUFFIX ".so") diff --git a/libsrc/csg/csg.hpp b/libsrc/csg/csg.hpp index 8ddf8d50..d524903f 100644 --- a/libsrc/csg/csg.hpp +++ b/libsrc/csg/csg.hpp @@ -24,9 +24,9 @@ #include "csgeom.hpp" #include "csgparser.hpp" - #include "triapprox.hpp" #include "algprim.hpp" +#include "splinesurface.hpp" #include "brick.hpp" #include "spline3d.hpp" #include "manifold.hpp" diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index de45e713..d0869662 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -338,6 +338,20 @@ namespace netgen FaceDescriptor & fd = mesh.GetFaceDescriptor(k); fd.SetBCName ( mesh.GetBCNamePtr ( fd.BCProperty() - 1 ) ); } + + int bbccnt = 0; + for (int k = 0; k < geom.GetNSurf(); k++){ + auto splinesurf = dynamic_cast (geom.GetSurface(k)); + if (splinesurf) + { + for( int i=0; i< splinesurf->GetNSplines(); i++) + { + string bcname = *splinesurf->GetBCName(i); + if(bcname != "default"){ + mesh.SetCD2Name(bbccnt,bcname); bbccnt++; } + } + } + } //!! diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 487ab002..dec21fb7 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -233,6 +233,41 @@ DLL_HEADER void ExportCSG() })) ; + bp::class_,shared_ptr>> ("SplineCurve3d") + .def ("AddPoint", FunctionPointer + ([] (SplineGeometry<3> & self, double x, double y, double z) + { + self.geompoints.Append (GeomPoint<3> (Point<3> (x,y,z))); + return self.geompoints.Size()-1; + })) + .def ("AddSegment", FunctionPointer + ([] (SplineGeometry<3> & self, int i1, int i2) + { + self.splines.Append (new LineSeg<3> (self.geompoints[i1], self.geompoints[i2])); + })) + .def ("AddSegment", FunctionPointer + ([] (SplineGeometry<3> & self, int i1, int i2, int i3) + { + self.splines.Append (new SplineSeg3<3> (self.geompoints[i1], self.geompoints[i2], self.geompoints[i3])); + })) + ; + + bp::class_,boost::noncopyable> ("SplineSurface") + .def("AddPoint", FunctionPointer + ([] (SplineSurface & self, double x, double y, double z, bool hpref) + { + self.AppendPoint(Point<3>(x,y,z),hpref); + return self.GetNP()-1; + }), + (bp::arg("self"),bp::arg("x"),bp::arg("y"),bp::arg("z"),bp::arg("hpref")=false)) + .def("AddSegment", FunctionPointer + ([] (SplineSurface & self, int i1, int i2, string bcname) + { + auto str = new string(bcname); + self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str); + }), + (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default")) + ; #if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100) bp::register_ptr_to_python>(); @@ -412,9 +447,19 @@ DLL_HEADER void ExportCSG() }), (bp::arg("self"), bp::arg("surface"), bp::arg("solid")) ) - + .def("AddSplineSurface", FunctionPointer + ([] (CSGeometry & self, shared_ptr surf) + { + auto planes = surf->CreatePlanes(); + auto spsol = make_shared(new Solid(&*surf)); + for(auto plane : (*planes)){ + spsol = make_shared(SPSolid::SECTION,spsol,make_shared(new Solid(plane))); + } + spsol->AddSurfaces(self); + int tlonr = self.SetTopLevelObject(spsol->GetSolid(), &*surf); + }), + (bp::arg("self"), bp::arg("SplineSurface"))) - .def("CloseSurfaces", FunctionPointer ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, bp::list aslices ) { diff --git a/libsrc/csg/splinesurface.cpp b/libsrc/csg/splinesurface.cpp new file mode 100644 index 00000000..92c116ed --- /dev/null +++ b/libsrc/csg/splinesurface.cpp @@ -0,0 +1,115 @@ + +#include + +namespace netgen +{ +SplineSurface :: SplineSurface() : OneSurfacePrimitive() +{ ; } + +SplineSurface :: ~SplineSurface() { ; } + +void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref) +{ + geompoints.Append(GeomPoint<3>(p,reffac)); + geompoints.Last().hpref = hpref; +} + + void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname) + { + splines.Append(spline); + bcnames.Append(bcname); + } + +double SplineSurface :: CalcFunctionValue (const Point<3> & point) const +{ + auto v1 = splines[0]->GetTangent(0); + auto v2 = splines.Last()->GetTangent(1); + auto p1 = splines[0]->GetPoint(0); + auto n = Cross(v1,v2); + n.Normalize(); + return n * (point-p1); +} + +void SplineSurface :: CalcGradient (const Point<3> & point, Vec<3> & grad) const +{ + auto v1 = splines[0]->GetTangent(0); + auto v2 = splines.Last()->GetTangent(1); + auto p1 = splines[0]->GetPoint(0); + grad = Cross(v1,v2); + grad.Normalize(); +} + + double SplineSurface :: HesseNorm() const +{ + return 0; +} + + Point<3> SplineSurface :: GetSurfacePoint() const +{ + return splines[0]->GetPoint(0); +} + + + void SplineSurface :: CalcSpecialPoints(Array> & pts) const + { + for(auto pt : geompoints) + { + pts.Append(Point<3>(pt)); + } + } + + Array* SplineSurface :: CreatePlanes() const + { + auto planes = new Array(); + auto sol = new Solid(new Plane(splines[0]->GetPoint(0),-(splines[0]->GetTangent(0)))); + for(auto spline : splines) + { + planes->Append(new Plane(spline->GetPoint(0),-spline->GetTangent(0))); + } + return planes; + } + + + void SplineSurface :: Print(ostream & str) const +{ + str << "SplineSurface " << endl; +} + +INSOLID_TYPE SplineSurface :: BoxInSolid(const BoxSphere<3> & box) const +{ + + int i; + double val; + + val = CalcFunctionValue (box.Center()); + if (val > box.Diam() / 2) return IS_OUTSIDE; + if (val < -box.Diam() / 2) return IS_INSIDE; + + Vec<3> n = GetNormalVector(box.Center()); + double cx = n[0]; + double cy = n[1]; + double cz = n[2]; + + if (val > 0) + { + Vec<3> vdiag = box.PMax() - box.PMin(); + double modify = (vdiag(0) * fabs (cx) + + vdiag(1) * fabs (cy) + + vdiag(2) * fabs (cz) ) / 2; + + if (val - modify < 0) + return DOES_INTERSECT; + return IS_OUTSIDE; + } + else + { + Vec<3> vdiag = box.PMax() - box.PMin(); + double modify = (vdiag(0) * fabs (cx) + + vdiag(1) * fabs (cy) + + vdiag(2) * fabs (cz) ) / 2; + if (val + modify > 0) + return DOES_INTERSECT; + return IS_INSIDE; + } +} +} diff --git a/libsrc/csg/splinesurface.cpp~ b/libsrc/csg/splinesurface.cpp~ new file mode 100644 index 00000000..1f6605dd --- /dev/null +++ b/libsrc/csg/splinesurface.cpp~ @@ -0,0 +1,92 @@ + +#include +#include +#include "splinesurface.hpp" + + +namespace netgen +{ +SplineSurface :: SplineSurface() : OneSurfacePrimitive() +{ ; } + +SplineSurface :: ~SplineSurface() { ; } + +void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref) +{ + geompoints.Append(GeomPoint<3>(p,reffac)); + geompoints.Last().hpref = hpref; +} + +double SplineSurface :: CalcFunctionValue (const Point<3> & point) const +{ + auto v1 = splines[0]->GetTangent(0); + auto v2 = splines.Last()->GetTangent(1); + auto p1 = splines[0]->GetPoint(0); + auto n = Cross(v1,v2); + n.Normalize(); + return n * (point-p1); +} + +void SplineSurface :: CalcGradient (const Point<3> & point, Vec<3> & grad) const +{ + auto v1 = splines[0]->GetTangent(0); + auto v2 = splines.Last()->GetTangent(1); + auto p1 = splines[0]->GetPoint(0); + grad = Cross(v1,v2); + grad.Normalize(); +} + + double SplineSurface :: HesseNorm() const +{ + return 0; +} + + Point<3> SplineSurface :: GetSurfacePoint() const +{ + return splines[0]->GetPoint(0); +} + + void SplineSurface :: Print(ostream & str) const +{ + str << "SplineSurface " << endl; +} + +INSOLID_TYPE SplineSurface :: BoxInSolid(const BoxSphere<3> & box) const +{ + +int i; +double val; + + + val = CalcFunctionValue (box.Center()); + if (val > box.Diam() / 2) return IS_OUTSIDE; + if (val < -box.Diam() / 2) return IS_INSIDE; + +Vec<3> n = GetNormalVector(box.Center()); +double cx = n[0]; +double cy = n[1]; +double cz = n[2]; + + if (val > 0) + { + Vec<3> vdiag = box.PMax() - box.PMin(); + double modify = (vdiag(0) * fabs (cx) + + vdiag(1) * fabs (cy) + + vdiag(2) * fabs (cz) ) / 2; + + if (val - modify < 0) + return DOES_INTERSECT; + return IS_OUTSIDE; + } + else + { + Vec<3> vdiag = box.PMax() - box.PMin(); + double modify = (vdiag(0) * fabs (cx) + + vdiag(1) * fabs (cy) + + vdiag(2) * fabs (cz) ) / 2; + if (val + modify > 0) + return DOES_INTERSECT; + return IS_INSIDE; +} +} +} diff --git a/libsrc/csg/splinesurface.hpp b/libsrc/csg/splinesurface.hpp new file mode 100644 index 00000000..63efa2c1 --- /dev/null +++ b/libsrc/csg/splinesurface.hpp @@ -0,0 +1,47 @@ +#ifndef FILE_SPLINESURFACE +#define FILE_SPLINESURFACE + + +namespace netgen +{ + class SplineSurface : public OneSurfacePrimitive + { + protected: + Array> geompoints; + Array*> splines; + Array bcnames; + + public: + SplineSurface(); + virtual ~SplineSurface(); + + const Array*> & GetSplines() const { return splines; } + int GetNSplines() const { return splines.Size(); } + string GetSplineType(const int i) const { return splines[i]->GetType(); } + SplineSeg<3> & GetSpline(const int i) { return *splines[i]; } + const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; } + int GetNP() const { return geompoints.Size(); } + const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; } + string* GetBCName(int i) const { return bcnames[i]; } + + DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false); + void AppendSegment(SplineSeg<3>* spline, string* bcname); + Array* CreatePlanes() const; + + virtual double CalcFunctionValue (const Point<3> & point) const; + virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const; + virtual double HesseNorm () const; + virtual Point<3> GetSurfacePoint () const; + virtual void CalcSpecialPoints(Array> & pts) const; + + virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const; + + + virtual void Print (ostream & str) const; + + + }; + +} + +#endif diff --git a/libsrc/csg/splinesurface.hpp~ b/libsrc/csg/splinesurface.hpp~ new file mode 100644 index 00000000..2f4127a0 --- /dev/null +++ b/libsrc/csg/splinesurface.hpp~ @@ -0,0 +1,46 @@ +#ifndef FILE_SPLINESURFACE +#define FILE_SPLINESURFACE + + +namespace netgen +{ + class SplineSurface : public OneSurfacePrimitive + { + protected: + Array> geompoints; + Array*> splines; + Array bcnames; + + public: + SplineSurface(); + virtual ~SplineSurface(); + + const Array*> & GetSplines() const { return splines; } + int GetNSplines() const { return splines.Size(); } + string GetSplineType(const int i) const { return splines[i]->GetType(); } + SplineSeg<3> & GetSpline(const int i) { return *splines[i]; } + const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; } + int GetNP() const { return geompoints.Size(); } + const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; } + + DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false); + void AppendSegment(SplineSeg<3>* spline, string* bcname); + Array* CreatePlanes() const; + + virtual double CalcFunctionValue (const Point<3> & point) const; + virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const; + virtual double HesseNorm () const; + virtual Point<3> GetSurfacePoint () const; + virtual void CalcSpecialPoints(Array> & pts) const; + + virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const; + + + virtual void Print (ostream & str) const; + + + }; + +} + +#endif diff --git a/libsrc/include/nginterface.h b/libsrc/include/nginterface.h index 1ffe2793..918f445d 100644 --- a/libsrc/include/nginterface.h +++ b/libsrc/include/nginterface.h @@ -125,6 +125,9 @@ extern "C" { DLL_HEADER char * Ng_GetBCNumBCName (int bcnr); //void Ng_GetBCNumBCName (int bcnr, char * name); + // Get BCName for bc-number of co dim 2 + DLL_HEADER char * Ng_GetCD2NumCD2Name (int cd2nr); + // Get normal vector of surface element node // DLL_HEADER void Ng_GetNormalVector (int sei, int locpi, double * nv); diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index 563e0b76..206e703d 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -542,6 +542,11 @@ char * Ng_GetBCNumBCName (int bcnr) return const_cast(mesh->GetBCName(bcnr).c_str()); } +char * Ng_GetCD2NumCD2Name (int cd2nr) +{ + return const_cast(mesh->GetCD2Name(cd2nr).c_str()); +} + // Inefficient (but maybe safer) version: //void Ng_GetBCNumBCName (int bcnr, char * name) diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index 92c1bbee..d3e8fc83 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -456,6 +456,22 @@ namespace netgen } } + template <> DLL_HEADER void Ngx_Mesh :: + ElementTransformation<1,3> (int elnr, + const double * xi, + double * x, + double * dxdxi) const + { + Point<3> xg; + Vec<3> dx; + mesh->GetCurvedElements().CalcSegmentTransformation(xi[0],elnr,xg,dx); + if(x) + for(int i=0;i<3;i++) x[i] = xg(i); + + if(dxdxi) + for(int i=0;i<3;i++) dxdxi[i] = dx(i); + } + template <> DLL_HEADER void Ngx_Mesh :: ElementTransformation<2,2> (int elnr, diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index f190e683..18fb1635 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -40,6 +40,7 @@ namespace netgen geomtype = NO_GEOM; bcnames.SetSize(0); + cd2names.SetSize(0); #ifdef PARALLEL paralleltop = new ParallelMeshTopology (*this); @@ -72,6 +73,9 @@ namespace netgen for (int i = 0; i < bcnames.Size(); i++ ) delete bcnames[i]; + for (int i = 0; i < cd2names.Size(); i++) + delete cd2names[i]; + #ifdef PARALLEL delete paralleltop; #endif @@ -94,6 +98,11 @@ namespace netgen if ( mesh2.bcnames[i] ) bcnames[i] = new string ( *mesh2.bcnames[i] ); else bcnames[i] = 0; + cd2names.SetSize(mesh2.cd2names.Size()); + for (int i=0; i < mesh2.cd2names.Size(); i++) + if (mesh2.cd2names[i]) cd2names[i] = new string(*mesh2.cd2names[i]); + else cd2names[i] = 0; + return *this; } @@ -126,6 +135,8 @@ namespace netgen for ( int i = 0; i < bcnames.Size(); i++ ) if ( bcnames[i] ) delete bcnames[i]; + for (int i= 0; i< cd2names.Size(); i++) + if (cd2names[i]) delete cd2names[i]; #ifdef PARALLEL delete paralleltop; @@ -610,6 +621,17 @@ namespace netgen outfile << i+1 << "\t" << GetBCName(i) << endl; outfile << endl << endl; } + int cntcd2names = 0; + for (int ii = 0; ii> n; + Array cd2nrs(n); + SetNCD2Names(n); + for( i=1; i<=n; i++) + { + string nextcd2name; + infile >> cd2nrs[i-1] >> nextcd2name; + cd2names[cd2nrs[i-1]-1] = new string(nextcd2name); + } + if (GetDimension() == 2) + { + throw NgException("co dim 2 elements not implemented for dimension 2"); + } + else + { + for (i = 1; i<= GetNSeg(); i++) + { + Segment & seg = LineSegment(i); + if ( seg.cd2i <= n ) + seg.SetBCName (cd2names[seg.cd2i-1]); + else + seg.SetBCName(0); + } + } + } + if (strcmp (str, "singular_points") == 0) { infile >> n; @@ -5714,6 +5762,46 @@ namespace netgen return defaultstring; } + void Mesh :: SetNCD2Names( int ncd2n ) + { + if (cd2names.Size()) + for(int i=0; i= cd2names.Size()) + { + int oldsize = cd2names.Size(); + cd2names.SetSize(cd2nr+1); + for(int i= oldsize; i<= cd2nr; i++) + cd2names[i] = nullptr; + } + if (cd2names[cd2nr]) delete cd2names[cd2nr]; + if (abcname != "default") + cd2names[cd2nr] = new string(abcname); + else + cd2names[cd2nr] = nullptr; + } + + const string & Mesh :: GetCD2Name (int cd2nr) const + { + static string defaultstring = "default"; + if (!cd2names.Size()) + return defaultstring; + + if (cd2nr < 0 || cd2nr >= cd2names.Size()) + throw NgException ("illegal bc-number"); + + if (cd2names[cd2nr]) + return *cd2names[cd2nr]; + else + return defaultstring; + } + void Mesh :: SetUserData(const char * id, Array & data) { if(userdata_int.Used(id)) diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index 3fce6709..a41b2b33 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -87,6 +87,9 @@ namespace netgen /// labels for boundary conditions Array bcnames; + /// labels for co dim 2 bboundary conditions + Array cd2names; + /// Periodic surface, close surface, etc. identifications Identifications * ident; @@ -588,6 +591,11 @@ namespace netgen const string & GetBCName ( int bcnr ) const; + DLL_HEADER void SetNCD2Names (int ncd2n); + DLL_HEADER void SetCD2Name (int cd2nr, const string & abcname); + + const string & GetCD2Name (int cd2nr ) const; + string * GetBCNamePtr (int bcnr) const { return bcnr < bcnames.Size() ? bcnames[bcnr] : nullptr; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 3687e240..78598a44 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -826,7 +826,9 @@ namespace netgen unsigned int seginfo:2; /// surface decoding index - int si; + int si; + /// co dim 2 deconding index + int cd2i; /// domain number inner side int domin; /// domain number outer side From 6134717796f85fd4da1c3667dd1b6d15c4519b75 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 5 Oct 2016 19:48:18 +0200 Subject: [PATCH 2/9] setting bboundary names --- libsrc/csg/edgeflw.cpp | 13 ++++++++++++- libsrc/csg/genmesh.cpp | 15 --------------- libsrc/csg/splinesurface.cpp | 15 +++++++++++++++ libsrc/csg/splinesurface.hpp | 1 + libsrc/include/nginterface_v2_impl.hpp | 10 ++++++++-- libsrc/meshing/meshclass.cpp | 5 +++-- libsrc/meshing/meshtype.cpp | 8 ++++++++ libsrc/meshing/meshtype.hpp | 6 ++++++ 8 files changed, 53 insertions(+), 20 deletions(-) diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index 5a5e22a6..b060d683 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -485,7 +485,18 @@ namespace netgen layer, mesh); } - + + (*testout) << "refedges size: " << refedges.Size() << endl; + for(int i=0; i(geometry.GetSurface(refedges[i].surfnr1)); + if(splinesurface) + { + auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); + mesh.SetCD2Name(refedges[i].edgenr-1,*name); + } + } /* // not available ... diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index d0869662..9c5e78bf 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -339,21 +339,6 @@ namespace netgen fd.SetBCName ( mesh.GetBCNamePtr ( fd.BCProperty() - 1 ) ); } - int bbccnt = 0; - for (int k = 0; k < geom.GetNSurf(); k++){ - auto splinesurf = dynamic_cast (geom.GetSurface(k)); - if (splinesurf) - { - for( int i=0; i< splinesurf->GetNSplines(); i++) - { - string bcname = *splinesurf->GetBCName(i); - if(bcname != "default"){ - mesh.SetCD2Name(bbccnt,bcname); bbccnt++; } - } - } - } - - //!! for (int k = 1; k <= mesh.GetNFD(); k++) diff --git a/libsrc/csg/splinesurface.cpp b/libsrc/csg/splinesurface.cpp index 92c116ed..46fd7a73 100644 --- a/libsrc/csg/splinesurface.cpp +++ b/libsrc/csg/splinesurface.cpp @@ -19,6 +19,21 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const splines.Append(spline); bcnames.Append(bcname); } + + string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const + { + (*testout) << "segment: " << p1 << ", " << p2 << endl; + for(int i=0; iGetPoint(0) << ", " << splines[i]->GetPoint(1) << endl; + if (((splines[i]->GetPoint(0)-p1).Length()<1e-12 && (splines[i]->GetPoint(1)-p2).Length() < 1e-12) || ((splines[i]->GetPoint(0)-p2).Length() < 1e-12 && (splines[i]->GetPoint(1)-p1).Length() < 1e-12)) + { + (*testout) << "return bcname: " << *bcnames[i] << endl; + return bcnames[i]; + } + } + return new string("default"); + } double SplineSurface :: CalcFunctionValue (const Point<3> & point) const { diff --git a/libsrc/csg/splinesurface.hpp b/libsrc/csg/splinesurface.hpp index 63efa2c1..70809163 100644 --- a/libsrc/csg/splinesurface.hpp +++ b/libsrc/csg/splinesurface.hpp @@ -23,6 +23,7 @@ namespace netgen int GetNP() const { return geompoints.Size(); } const GeomPoint<3> & GetPoint(int i) const { return geompoints[i]; } string* GetBCName(int i) const { return bcnames[i]; } + string* GetBCNameOf(Point<3> p1, Point<3> p2) const; DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false); void AppendSegment(SplineSeg<3>* spline, string* bcname); diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index b7689b8a..dd94d02f 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -13,7 +13,10 @@ NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<0> (int nr) const template <> NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const { - return (*mesh)[SegmentIndex(nr)].si; + if(mesh->GetDimension()==3) + return (*mesh)[SegmentIndex(nr)].cd2i; + else + return (*mesh)[SegmentIndex(nr)].si; } template <> @@ -63,7 +66,10 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const Ng_Element ret; ret.type = NG_ELEMENT_TYPE(el.GetType()); - ret.index = el.si; + if(mesh->GetDimension()==2) + ret.index = el.si; + else + ret.index = el.edgenr; ret.points.num = el.GetNP(); ret.points.ptr = (int*)&(el[0]); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 18fb1635..ba19a36c 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5773,6 +5773,7 @@ namespace netgen void Mesh :: SetCD2Name ( int cd2nr, const string & abcname ) { + (*testout) << "setCD2Name on edge " << cd2nr << " to " << abcname << endl; if (cd2nr >= cd2names.Size()) { int oldsize = cd2names.Size(); @@ -5780,7 +5781,7 @@ namespace netgen for(int i= oldsize; i<= cd2nr; i++) cd2names[i] = nullptr; } - if (cd2names[cd2nr]) delete cd2names[cd2nr]; + //if (cd2names[cd2nr]) delete cd2names[cd2nr]; if (abcname != "default") cd2names[cd2nr] = new string(abcname); else @@ -5794,7 +5795,7 @@ namespace netgen return defaultstring; if (cd2nr < 0 || cd2nr >= cd2names.Size()) - throw NgException ("illegal bc-number"); + return defaultstring; if (cd2names[cd2nr]) return *cd2names[cd2nr]; diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index 2dadf981..acad5a7f 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2418,6 +2418,14 @@ namespace netgen return s; } + string EdgeDescriptor :: default_bcname = "default"; + void EdgeDescriptor :: SetBCName (string * bcn) + { + if(bcn) + bcname = bcn; + else + bcn = &default_bcname; + } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 78598a44..e4963934 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1005,6 +1005,9 @@ namespace netgen { int tlosurf; int surfnr[2]; + int bcprop; + static string default_bcname; + string* bcname = &default_bcname; public: EdgeDescriptor () : tlosurf(-1) @@ -1015,6 +1018,9 @@ namespace netgen int TLOSurface() const { return tlosurf; } void SetTLOSurface (int nr) { tlosurf = nr; } + int BCProperty() const { return bcprop; } + void SetBCProperty (int bc) { bcprop = bc; } + void SetBCName (string* bcn); }; From b8bf194fcfe300e584d36ad5a367ac272269ea51 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 11 Oct 2016 14:10:36 +0200 Subject: [PATCH 3/9] some changes --- libsrc/csg/algprim.hpp | 1 + libsrc/csg/edgeflw.cpp | 23 ++++++++++++-- libsrc/csg/genmesh.cpp | 6 ++++ libsrc/csg/splinesurface.cpp | 43 ++++++++++++++++---------- libsrc/csg/splinesurface.hpp | 6 ++-- libsrc/include/nginterface_v2_impl.hpp | 2 +- libsrc/meshing/meshclass.cpp | 1 + libsrc/meshing/meshfunc.cpp | 2 +- libsrc/meshing/meshtype.cpp | 12 ------- libsrc/meshing/meshtype.hpp | 6 ---- 10 files changed, 60 insertions(+), 42 deletions(-) diff --git a/libsrc/csg/algprim.hpp b/libsrc/csg/algprim.hpp index 0459108a..9ca93cdd 100644 --- a/libsrc/csg/algprim.hpp +++ b/libsrc/csg/algprim.hpp @@ -53,6 +53,7 @@ namespace netgen /// A Plane (i.e., the plane and everything behind it). class Plane : public QuadraticSurface { + protected: /// a point in the plane Point<3> p; /// outward normal vector diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index b060d683..4541567c 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -486,15 +486,13 @@ namespace netgen mesh); } - (*testout) << "refedges size: " << refedges.Size() << endl; for(int i=0; i(geometry.GetSurface(refedges[i].surfnr1)); if(splinesurface) { auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); - mesh.SetCD2Name(refedges[i].edgenr-1,*name); + mesh.SetCD2Name(refedges[i].edgenr,*name); } } @@ -1168,7 +1166,23 @@ namespace netgen refedges.Elem(hi).domout = i; } else + { refedges.Elem(hi).tlosurf = i; + for(int kk = 0; kk < geometry.GetNTopLevelObjects(); kk++) + { + auto othersolid = geometry.GetTopLevelObject(kk)->GetSolid(); + auto othersurf = geometry.GetTopLevelObject(kk)->GetSurface(); + if(!othersurf) + { + if(othersolid->IsIn(edgepoints[0]) && + othersolid->IsIn(edgepoints[edgepoints.Size()-1])) + { + refedges.Elem(hi).domin = kk; + refedges.Elem(hi).domout = kk; + } + } + } + } if(pre_ok[k-1]) edges_priority[hi-1] = 1; @@ -1386,6 +1400,7 @@ namespace netgen seg.surfnr2 = refedges.Get(k).surfnr2; seg.seginfo = 0; if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; + *testout << "add segment at 3" << endl; mesh.AddSegment (seg); //(*testout) << "add seg " << mesh[seg.p1] << "-" << mesh[seg.p2] << endl; //(*testout) << "refedge " << k << " surf1 " << seg.surfnr1 << " surf2 " << seg.surfnr2 << " inv " << refedgesinv.Get(k) << endl; @@ -1554,6 +1569,7 @@ namespace netgen seg.surfnr2 = refedges.Get(k).surfnr2; seg.seginfo = 0; if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; + *testout << "add segment at 1" << endl; mesh.AddSegment (seg); // (*testout) << "add seg " << seg[0] << "-" << seg[1] << endl; } @@ -1685,6 +1701,7 @@ namespace netgen seg.surfnr2 = refedges.Get(k).surfnr2; seg.seginfo = 0; if (k == 1) seg.seginfo = refedgesinv.Get(k) ? 2 : 1; + *testout << "add segment at 2" << endl; mesh.AddSegment (seg); // (*testout) << "copy seg " << seg[0] << "-" << seg[1] << endl; #ifdef DEVELOP diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index 9c5e78bf..07fa41d2 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -98,6 +98,12 @@ namespace netgen //(*testout) << " to " << mesh.LineSegment(i).si << endl; } + for(int k = 1; k<=mesh.GetNFD(); k++) + { + *testout << "face: " << k << endl + << "FD: " << mesh.GetFaceDescriptor(k) << endl; + } + if (geom.identifications.Size()) { PrintMessage (3, "Find Identifications"); diff --git a/libsrc/csg/splinesurface.cpp b/libsrc/csg/splinesurface.cpp index 46fd7a73..1009f41f 100644 --- a/libsrc/csg/splinesurface.cpp +++ b/libsrc/csg/splinesurface.cpp @@ -3,7 +3,7 @@ namespace netgen { -SplineSurface :: SplineSurface() : OneSurfacePrimitive() + SplineSurface :: SplineSurface() : OneSurfacePrimitive() { ; } SplineSurface :: ~SplineSurface() { ; } @@ -34,7 +34,30 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const } return new string("default"); } - + + Array* SplineSurface :: CreatePlanes() const + { + auto planes = new Array(); + auto sol = new Solid(new Plane(splines[0]->GetPoint(0),Cross(splines[0]->GetTangent(0),GetNormalVector(splines[0]->GetPoint(0))))); + for(auto spline : splines) + { + planes->Append(new Plane(spline->GetPoint(0),-spline->GetTangent(0))); + } + return planes; + } + + void SplineSurface :: Print(ostream & str) const +{ + str << "SplineSurface " << endl; +} + + void SplineSurface :: Project(Point<3> & p3d) const + { + double val = CalcFunctionValue(p3d); + p3d -= val * GetNormalVector(p3d); + } + + double SplineSurface :: CalcFunctionValue (const Point<3> & point) const { auto v1 = splines[0]->GetTangent(0); @@ -73,23 +96,8 @@ void SplineSurface :: CalcGradient (const Point<3> & point, Vec<3> & grad) const } } - Array* SplineSurface :: CreatePlanes() const - { - auto planes = new Array(); - auto sol = new Solid(new Plane(splines[0]->GetPoint(0),-(splines[0]->GetTangent(0)))); - for(auto spline : splines) - { - planes->Append(new Plane(spline->GetPoint(0),-spline->GetTangent(0))); - } - return planes; - } - void SplineSurface :: Print(ostream & str) const -{ - str << "SplineSurface " << endl; -} - INSOLID_TYPE SplineSurface :: BoxInSolid(const BoxSphere<3> & box) const { @@ -127,4 +135,5 @@ INSOLID_TYPE SplineSurface :: BoxInSolid(const BoxSphere<3> & box) const return IS_INSIDE; } } + } diff --git a/libsrc/csg/splinesurface.hpp b/libsrc/csg/splinesurface.hpp index 70809163..9211049f 100644 --- a/libsrc/csg/splinesurface.hpp +++ b/libsrc/csg/splinesurface.hpp @@ -4,7 +4,7 @@ namespace netgen { - class SplineSurface : public OneSurfacePrimitive + class SplineSurface : public OneSurfacePrimitive { protected: Array> geompoints; @@ -28,7 +28,8 @@ namespace netgen DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false); void AppendSegment(SplineSeg<3>* spline, string* bcname); Array* CreatePlanes() const; - + + virtual void Project (Point<3> & p3d) const; virtual double CalcFunctionValue (const Point<3> & point) const; virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const; virtual double HesseNorm () const; @@ -36,6 +37,7 @@ namespace netgen virtual void CalcSpecialPoints(Array> & pts) const; virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const; + virtual void Print (ostream & str) const; diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index dd94d02f..e6918e7e 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -14,7 +14,7 @@ template <> NGX_INLINE DLL_HEADER int Ngx_Mesh :: GetElementIndex<1> (int nr) const { if(mesh->GetDimension()==3) - return (*mesh)[SegmentIndex(nr)].cd2i; + return (*mesh)[SegmentIndex(nr)].edgenr; else return (*mesh)[SegmentIndex(nr)].si; } diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index ba19a36c..36d07b1d 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5773,6 +5773,7 @@ namespace netgen void Mesh :: SetCD2Name ( int cd2nr, const string & abcname ) { + cd2nr++; (*testout) << "setCD2Name on edge " << cd2nr << " to " << abcname << endl; if (cd2nr >= cd2names.Size()) { diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index 8440cc1f..54f06a65 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -86,7 +86,7 @@ namespace netgen { const Element2d & el = mesh3d[sei]; if (el.IsDeleted() ) continue; - + if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k || mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k) diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index acad5a7f..b4476e00 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2418,18 +2418,6 @@ namespace netgen return s; } - string EdgeDescriptor :: default_bcname = "default"; - void EdgeDescriptor :: SetBCName (string * bcn) - { - if(bcn) - bcname = bcn; - else - bcn = &default_bcname; - } - - - - Identifications :: Identifications (Mesh & amesh) : mesh(amesh) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index e4963934..78598a44 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1005,9 +1005,6 @@ namespace netgen { int tlosurf; int surfnr[2]; - int bcprop; - static string default_bcname; - string* bcname = &default_bcname; public: EdgeDescriptor () : tlosurf(-1) @@ -1018,9 +1015,6 @@ namespace netgen int TLOSurface() const { return tlosurf; } void SetTLOSurface (int nr) { tlosurf = nr; } - int BCProperty() const { return bcprop; } - void SetBCProperty (int bc) { bcprop = bc; } - void SetBCName (string* bcn); }; From ccde47d2c07e5b5c9a6b5fb38fc38e407d1938d7 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Sun, 16 Oct 2016 09:45:16 +0200 Subject: [PATCH 4/9] added multielementtransformation --- libsrc/interface/nginterface_v2.cpp | 63 +++++++++++++++++++++++++++++ libsrc/meshing/curvedelems.cpp | 1 + 2 files changed, 64 insertions(+) diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index d3e8fc83..d7548127 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -535,6 +535,17 @@ namespace netgen if (dxdxi) dxdxi[0] = dx(0); } + template <> DLL_HEADER void Ngx_Mesh :: + ElementTransformation<0,2> (int elnr, + const double *xi, + double * x, + double * dxdxi) const + { + PointIndex pnum = mesh->pointelements[elnr].pnum; + if (x) + for (int i = 0; i< 2; i++) x[i] = (*mesh)[pnum](i); + } + template <> DLL_HEADER void Ngx_Mesh :: ElementTransformation<0,1> (int elnr, @@ -581,6 +592,15 @@ namespace netgen mesh->GetCurvedElements().CalcMultiPointSurfaceTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi); } + template <> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<1,3> (int elnr, int npts, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi) const + { + mesh->GetCurvedElements().CalcMultiPointSegmentTransformation<3> (elnr, npts, xi, sxi, x, sx, dxdxi, sdxdxi); + } + template <> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<1,2> (int elnr, int npts, const double * xi, size_t sxi, @@ -600,6 +620,16 @@ namespace netgen ElementTransformation<1,1> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi); } + template <> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<0,2> (int elnr, int npts, + const double * xi, size_t sxi, + double * x, size_t sx, + double * dxdxi, size_t sdxdxi) const + { + for (int i = 0; i < npts; i++) + ElementTransformation<0,2> (elnr, xi + i*sxi, x+i*sx, dxdxi+i*sdxdxi); + } + template <> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<0,1> (int elnr, int npts, @@ -698,6 +728,15 @@ namespace netgen */ } + template<> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<0,2> (int elnr, int npts, + const __m256d *xi, size_t sxi, + __m256d * x, size_t sx, + __m256d * dxdxi, size_t sdxdxi) const + { + cout << "MultiElementtransformation<0,2> simd not implemented" << endl; + } + template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<0,1> (int elnr, int npts, const __m256d * xi, size_t sxi, @@ -707,6 +746,30 @@ namespace netgen cout << "multi-eltrafo simd called, 0,1,simd" << endl; } + template<> DLL_HEADER void Ngx_Mesh :: + MultiElementTransformation<1,3> (int elnr, int npts, + const __m256d * xi, size_t sxi, + __m256d * x, size_t sx, + __m256d * dxdxi, size_t sdxdxi) const + { + double hxi[4][1]; + double hx[4][3]; + double hdxdxi[4][3]; + for (int j = 0; j<4;j++) + hxi[j][0] = ((double*)&(xi[0]))[j]; + MultiElementTransformation<1,3> (elnr, 4, &hxi[0][0], 1, &hx[0][0], 2, &hdxdxi[0][0],4); + for(int j=0; j<4; j++) + for(int k=0; k<3; k++) + ((double*)&(x[k]))[j] = hx[j][k]; + for(int j=0; j< 4; j++) + for (int k = 0; k<3; k++) + ((double*) & (dxdxi[k]))[j] = hdxdxi[j][k]; + + xi += sxi; + x += sx; + dxdxi += sdxdxi; + } + template<> DLL_HEADER void Ngx_Mesh :: MultiElementTransformation<1,2> (int elnr, int npts, const __m256d * xi, size_t sxi, diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 6f05ed86..14b0d4df 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -3686,6 +3686,7 @@ namespace netgen } } + template void CurvedElements :: CalcMultiPointSegmentTransformation<2> (SegmentIndex elnr, int npts, const double * xi, size_t sxi, From a6ea18d07d93226fe77dbdd54631cb14a96af14b Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Mon, 17 Oct 2016 17:31:09 +0200 Subject: [PATCH 5/9] fix bboundary condition numbering --- libsrc/meshing/meshclass.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 36d07b1d..b890399b 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -5773,7 +5773,7 @@ namespace netgen void Mesh :: SetCD2Name ( int cd2nr, const string & abcname ) { - cd2nr++; + cd2nr--; (*testout) << "setCD2Name on edge " << cd2nr << " to " << abcname << endl; if (cd2nr >= cd2names.Size()) { From 7a6de7b1dcc294af707256783dfbb661dbd6c62d Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 27 Oct 2016 15:41:08 +0200 Subject: [PATCH 6/9] further fixes --- libsrc/csg/edgeflw.cpp | 14 +++- libsrc/csg/python_csg.cpp | 37 ++++++++-- libsrc/csg/splinesurface.cpp | 136 +++++++++-------------------------- libsrc/csg/splinesurface.hpp | 41 +++++++++-- 4 files changed, 109 insertions(+), 119 deletions(-) diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index 4541567c..9ba3af9f 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -44,6 +44,7 @@ namespace netgen for (PointIndex pi : mesh.Points().Range()) meshpoint_tree->Insert (mesh[pi], pi); + // add all special points before edge points (important for periodic identification) // JS, Jan 2007 const double di=1e-7*geometry.MaxSize(); @@ -454,7 +455,6 @@ namespace netgen refedges[i].surfnr2 = geometry.GetSurfaceClassRepresentant(refedges[i].surfnr2); } - /* for (int i = oldnseg+1; i <= mesh.GetNSeg(); i++) for (int j = 1; j <= oldnseg; j++) @@ -485,7 +485,6 @@ namespace netgen layer, mesh); } - for(int i=0; i(geometry.GetSurface(refedges[i].surfnr1)); @@ -494,6 +493,17 @@ namespace netgen auto name = splinesurface->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); mesh.SetCD2Name(refedges[i].edgenr,*name); } + else + { + auto splinesurface2 = dynamic_cast(geometry.GetSurface(refedges[i].surfnr2)); + if(splinesurface2) + { + auto name = splinesurface2->GetBCNameOf(specpoints[startpoints.Get(refedges[i].edgenr)].p,specpoints[endpoints.Get(refedges[i].edgenr)].p); + mesh.SetCD2Name(refedges[i].edgenr,*name); + } + + } + } /* diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index dec21fb7..4ef7cd26 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -252,7 +252,28 @@ DLL_HEADER void ExportCSG() })) ; - bp::class_,boost::noncopyable> ("SplineSurface") + bp::class_,boost::noncopyable> ("SplineSurface", + "A surface for co dim 2 integrals on the splines", bp::no_init) + .def("__init__", bp::make_constructor (FunctionPointer + ([](const shared_ptr base, bp::list cuts) + { + auto primitive = dynamic_cast (base->GetSolid()->GetPrimitive()); + auto acuts = new Array(); + for(int i = 0; i> sps(cuts[i]); + if(!sps.check()) + throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!"); + auto sp = dynamic_cast(sps()->GetSolid()->GetPrimitive()); + if(sp) + acuts->Append(sp); + else + throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!"); + } + if(!primitive) + throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!"); + return make_shared(primitive,acuts); + }),bp::default_call_policies(),(bp::arg("base"), bp::arg("cuts")=bp::list()))) .def("AddPoint", FunctionPointer ([] (SplineSurface & self, double x, double y, double z, bool hpref) { @@ -261,12 +282,12 @@ DLL_HEADER void ExportCSG() }), (bp::arg("self"),bp::arg("x"),bp::arg("y"),bp::arg("z"),bp::arg("hpref")=false)) .def("AddSegment", FunctionPointer - ([] (SplineSurface & self, int i1, int i2, string bcname) + ([] (SplineSurface & self, int i1, int i2, string bcname, double maxh) { auto str = new string(bcname); - self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str); + self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str,maxh); }), - (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default")) + (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default", bp::arg("maxh")=-1.)) ; #if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100) @@ -450,13 +471,15 @@ DLL_HEADER void ExportCSG() .def("AddSplineSurface", FunctionPointer ([] (CSGeometry & self, shared_ptr surf) { - auto planes = surf->CreatePlanes(); + auto cuttings = surf->CreateCuttingSurfaces(); auto spsol = make_shared(new Solid(&*surf)); - for(auto plane : (*planes)){ - spsol = make_shared(SPSolid::SECTION,spsol,make_shared(new Solid(plane))); + for(auto cut : (*cuttings)){ + spsol = make_shared(SPSolid::SECTION,spsol,make_shared(new Solid(cut))); } spsol->AddSurfaces(self); int tlonr = self.SetTopLevelObject(spsol->GetSolid(), &*surf); + for(auto p : surf->GetPoints()) + self.AddUserPoint(p); }), (bp::arg("self"), bp::arg("SplineSurface"))) diff --git a/libsrc/csg/splinesurface.cpp b/libsrc/csg/splinesurface.cpp index 1009f41f..d66c8bff 100644 --- a/libsrc/csg/splinesurface.cpp +++ b/libsrc/csg/splinesurface.cpp @@ -3,47 +3,60 @@ namespace netgen { - SplineSurface :: SplineSurface() : OneSurfacePrimitive() -{ ; } - -SplineSurface :: ~SplineSurface() { ; } - void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref) { - geompoints.Append(GeomPoint<3>(p,reffac)); + auto pp = Point<3>(p); + geompoints.Append(GeomPoint<3>(pp,reffac)); geompoints.Last().hpref = hpref; } - void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname) + void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh) { - splines.Append(spline); - bcnames.Append(bcname); + splines.Append(spline); + bcnames.Append(bcname); + maxh.Append(amaxh); } string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const { - (*testout) << "segment: " << p1 << ", " << p2 << endl; + + double eps = 1e-5; for(int i=0; iGetPoint(0) << ", " << splines[i]->GetPoint(1) << endl; - if (((splines[i]->GetPoint(0)-p1).Length()<1e-12 && (splines[i]->GetPoint(1)-p2).Length() < 1e-12) || ((splines[i]->GetPoint(0)-p2).Length() < 1e-12 && (splines[i]->GetPoint(1)-p1).Length() < 1e-12)) + auto pp1 = Point<3>(splines[i]->GetPoint(0)); + Project(pp1); + auto pp2 = Point<3>(splines[i]->GetPoint(1)); + Project(pp2); + if (((pp1-p1).Length()* SplineSurface :: CreatePlanes() const + Array* SplineSurface :: CreateCuttingSurfaces() const { - auto planes = new Array(); - auto sol = new Solid(new Plane(splines[0]->GetPoint(0),Cross(splines[0]->GetTangent(0),GetNormalVector(splines[0]->GetPoint(0))))); - for(auto spline : splines) + auto cuttings = new Array(); + for (auto cut : *cuts) + cuttings->Append(cut); + for(int i = 0; iAppend(new Plane(spline->GetPoint(0),-spline->GetTangent(0))); + auto spline = splines[i]; + auto lineseg = dynamic_cast*>(spline); + auto p1 = Point<3>(spline->GetPoint(0)); + Project(p1); + auto p2 = Point<3>(spline->GetPoint(1)); + Project(p2); + auto vec = Vec<3>(p2)-Vec<3>(p1); + auto plane = new Plane(p1,-Cross(vec,baseprimitive->GetNormalVector(p1))); + if(maxh[i]>0) + { + plane->SetMaxH(maxh[i]); + } + cuttings->Append(plane); } - return planes; + return cuttings; } void SplineSurface :: Print(ostream & str) const @@ -51,89 +64,4 @@ void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const str << "SplineSurface " << endl; } - void SplineSurface :: Project(Point<3> & p3d) const - { - double val = CalcFunctionValue(p3d); - p3d -= val * GetNormalVector(p3d); - } - - -double SplineSurface :: CalcFunctionValue (const Point<3> & point) const -{ - auto v1 = splines[0]->GetTangent(0); - auto v2 = splines.Last()->GetTangent(1); - auto p1 = splines[0]->GetPoint(0); - auto n = Cross(v1,v2); - n.Normalize(); - return n * (point-p1); -} - -void SplineSurface :: CalcGradient (const Point<3> & point, Vec<3> & grad) const -{ - auto v1 = splines[0]->GetTangent(0); - auto v2 = splines.Last()->GetTangent(1); - auto p1 = splines[0]->GetPoint(0); - grad = Cross(v1,v2); - grad.Normalize(); -} - - double SplineSurface :: HesseNorm() const -{ - return 0; -} - - Point<3> SplineSurface :: GetSurfacePoint() const -{ - return splines[0]->GetPoint(0); -} - - - void SplineSurface :: CalcSpecialPoints(Array> & pts) const - { - for(auto pt : geompoints) - { - pts.Append(Point<3>(pt)); - } - } - - - -INSOLID_TYPE SplineSurface :: BoxInSolid(const BoxSphere<3> & box) const -{ - - int i; - double val; - - val = CalcFunctionValue (box.Center()); - if (val > box.Diam() / 2) return IS_OUTSIDE; - if (val < -box.Diam() / 2) return IS_INSIDE; - - Vec<3> n = GetNormalVector(box.Center()); - double cx = n[0]; - double cy = n[1]; - double cz = n[2]; - - if (val > 0) - { - Vec<3> vdiag = box.PMax() - box.PMin(); - double modify = (vdiag(0) * fabs (cx) + - vdiag(1) * fabs (cy) + - vdiag(2) * fabs (cz) ) / 2; - - if (val - modify < 0) - return DOES_INTERSECT; - return IS_OUTSIDE; - } - else - { - Vec<3> vdiag = box.PMax() - box.PMin(); - double modify = (vdiag(0) * fabs (cx) + - vdiag(1) * fabs (cy) + - vdiag(2) * fabs (cz) ) / 2; - if (val + modify > 0) - return DOES_INTERSECT; - return IS_INSIDE; - } -} - } diff --git a/libsrc/csg/splinesurface.hpp b/libsrc/csg/splinesurface.hpp index 9211049f..1f68ae73 100644 --- a/libsrc/csg/splinesurface.hpp +++ b/libsrc/csg/splinesurface.hpp @@ -4,19 +4,25 @@ namespace netgen { - class SplineSurface : public OneSurfacePrimitive + class SplineSurface : public OneSurfacePrimitive { protected: Array> geompoints; Array*> splines; Array bcnames; + Array maxh; + OneSurfacePrimitive* baseprimitive; + Array* cuts; public: - SplineSurface(); - virtual ~SplineSurface(); + SplineSurface(OneSurfacePrimitive* abaseprimitive, Array* acuts) : + OneSurfacePrimitive(), baseprimitive(abaseprimitive), cuts(acuts) + { ; } + virtual ~SplineSurface() { ; } const Array*> & GetSplines() const { return splines; } int GetNSplines() const { return splines.Size(); } + const Array>& GetPoints() const { return geompoints; } string GetSplineType(const int i) const { return splines[i]->GetType(); } SplineSeg<3> & GetSpline(const int i) { return *splines[i]; } const SplineSeg<3> & GetSpline(const int i) const { return *splines[i]; } @@ -26,9 +32,30 @@ namespace netgen string* GetBCNameOf(Point<3> p1, Point<3> p2) const; DLL_HEADER void AppendPoint(const Point<3> & p, const double reffac = 1., const bool hpref=false); - void AppendSegment(SplineSeg<3>* spline, string* bcname); - Array* CreatePlanes() const; + void AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh = -1); + OneSurfacePrimitive* GetBase() const { return baseprimitive; } + + Array* CreateCuttingSurfaces() const; + + + virtual void Project (Point<3> & p3d) const { baseprimitive->Project(p3d); } + virtual double CalcFunctionValue (const Point<3> & point) const + { return baseprimitive->CalcFunctionValue (point); } + virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const + { baseprimitive->CalcGradient (point,grad); } + virtual double HesseNorm () const + { return baseprimitive->HesseNorm(); } + virtual Point<3> GetSurfacePoint () const + { return baseprimitive->GetSurfacePoint(); } + virtual void CalcSpecialPoints(Array> & pts) const + { baseprimitive->CalcSpecialPoints(pts); } + + virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const + { return baseprimitive->BoxInSolid(box); } + + + /* virtual void Project (Point<3> & p3d) const; virtual double CalcFunctionValue (const Point<3> & point) const; virtual void CalcGradient (const Point<3> & point, Vec<3> & grad) const; @@ -38,13 +65,15 @@ namespace netgen virtual INSOLID_TYPE BoxInSolid(const BoxSphere<3> & box) const; - + */ virtual void Print (ostream & str) const; }; + + } #endif From 767b819e5d897a4cd1eae893dd747a0ef7d64f7f Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Fri, 28 Oct 2016 16:49:50 +0200 Subject: [PATCH 7/9] python cd2 functionality, mesh scaling --- libsrc/csg/edgeflw.cpp | 3 --- libsrc/meshing/meshclass.cpp | 4 +++- libsrc/meshing/meshclass.hpp | 1 + libsrc/meshing/meshtype.hpp | 2 ++ libsrc/meshing/python_mesh.cpp | 16 ++++++++++++++++ 5 files changed, 22 insertions(+), 4 deletions(-) diff --git a/libsrc/csg/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index 9ba3af9f..6afa4957 100644 --- a/libsrc/csg/edgeflw.cpp +++ b/libsrc/csg/edgeflw.cpp @@ -1410,7 +1410,6 @@ namespace netgen seg.surfnr2 = refedges.Get(k).surfnr2; seg.seginfo = 0; if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; - *testout << "add segment at 3" << endl; mesh.AddSegment (seg); //(*testout) << "add seg " << mesh[seg.p1] << "-" << mesh[seg.p2] << endl; //(*testout) << "refedge " << k << " surf1 " << seg.surfnr1 << " surf2 " << seg.surfnr2 << " inv " << refedgesinv.Get(k) << endl; @@ -1579,7 +1578,6 @@ namespace netgen seg.surfnr2 = refedges.Get(k).surfnr2; seg.seginfo = 0; if (k == 1) seg.seginfo = (refedgesinv.Get(k)) ? 2 : 1; - *testout << "add segment at 1" << endl; mesh.AddSegment (seg); // (*testout) << "add seg " << seg[0] << "-" << seg[1] << endl; } @@ -1711,7 +1709,6 @@ namespace netgen seg.surfnr2 = refedges.Get(k).surfnr2; seg.seginfo = 0; if (k == 1) seg.seginfo = refedgesinv.Get(k) ? 2 : 1; - *testout << "add segment at 2" << endl; mesh.AddSegment (seg); // (*testout) << "copy seg " << seg[0] << "-" << seg[1] << endl; #ifdef DEVELOP diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index b890399b..2423c03e 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -1120,7 +1120,7 @@ namespace netgen { Segment & seg = LineSegment(i); if ( seg.cd2i <= n ) - seg.SetBCName (cd2names[seg.cd2i-1]); + seg.SetBCName (cd2names[seg.edgenr-1]); else seg.SetBCName(0); } @@ -1407,7 +1407,9 @@ namespace netgen if(seg.surfnr2 >= 0) seg.surfnr2 = seg.surfnr2 + max_surfnr; seg[0] = seg[0] +oldnp; seg[1] = seg[1] +oldnp; + *testout << "old edgenr: " << seg.edgenr << endl; seg.edgenr = seg.edgenr + oldne; + *testout << "new edgenr: " << seg.edgenr << endl; seg.epgeominfo[1].edgenr = seg.epgeominfo[1].edgenr + oldne; AddSegment (seg); diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index a41b2b33..d09b7c98 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -595,6 +595,7 @@ namespace netgen DLL_HEADER void SetCD2Name (int cd2nr, const string & abcname); const string & GetCD2Name (int cd2nr ) const; + size_t GetNCD2Names() const { return cd2names.Size(); } string * GetBCNamePtr (int bcnr) const { return bcnr < bcnames.Size() ? bcnames[bcnr] : nullptr; } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 78598a44..80be7fa1 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -251,6 +251,8 @@ namespace netgen singular = 0; } + void Scale(double factor) { *testout << "before: " << x[0] << endl; x[0] *= factor; x[1] *= factor; x[2] *= factor; *testout << "after: " << x[0] << endl;} + int GetLayer() const { return layer; } POINTTYPE Type() const { return type; } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index eb14c78d..fd6a8c37 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -216,6 +216,7 @@ DLL_HEADER void ExportNetgenMeshing() tmp->surfnr1 = bp::extract(surfaces[0]); tmp->surfnr2 = bp::extract(surfaces[1]); } + tmp->edgenr = index; return tmp; }), bp::default_call_policies(), @@ -242,6 +243,10 @@ DLL_HEADER void ExportNetgenMeshing() li.append (self.surfnr2); return li; })) + .add_property("index", FunctionPointer([](const Segment &self) -> size_t + { + return self.edgenr; + })) ; @@ -465,6 +470,11 @@ DLL_HEADER void ExportNetgenMeshing() .def ("SetMaterial", &Mesh::SetMaterial) .def ("GetMaterial", FunctionPointer([](Mesh & self, int domnr) { return string(self.GetMaterial(domnr)); })) + + .def ("SetCD2Name", &Mesh::SetCD2Name) + .def ("GetCD2Name", FunctionPointer([](Mesh & self, int nr) -> string + { return self.GetCD2Name(nr); })) + .def ("GetNCD2Names", &Mesh::GetNCD2Names) .def ("GenerateVolumeMesh", FunctionPointer ([](Mesh & self) @@ -564,6 +574,12 @@ DLL_HEADER void ExportNetgenMeshing() GenerateBoundaryLayer (self, blp); } )) + + .def ("Scale", FunctionPointer([](Mesh & self, double factor) + { + for(auto i = 0; i Date: Sun, 6 Nov 2016 18:25:38 +0100 Subject: [PATCH 8/9] pyb11 fixes --- libsrc/csg/python_csg.cpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 8fc33039..9cd9390d 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -254,7 +254,7 @@ DLL_HEADER void ExportCSG(py::module &m) })) ; - bp::class_,shared_ptr>> ("SplineCurve3d") + py::class_,shared_ptr>> (m,"SplineCurve3d") .def ("AddPoint", FunctionPointer ([] (SplineGeometry<3> & self, double x, double y, double z) { @@ -273,16 +273,15 @@ DLL_HEADER void ExportCSG(py::module &m) })) ; - bp::class_,boost::noncopyable> ("SplineSurface", - "A surface for co dim 2 integrals on the splines", bp::no_init) - .def("__init__", bp::make_constructor (FunctionPointer - ([](const shared_ptr base, bp::list cuts) + py::class_> (m, "SplineSurface", + "A surface for co dim 2 integrals on the splines") + .def("__init__", FunctionPointer ([](const shared_ptr base, py::list cuts) { auto primitive = dynamic_cast (base->GetSolid()->GetPrimitive()); auto acuts = new Array(); - for(int i = 0; i> sps(cuts[i]); + py::extract> sps(cuts[i]); if(!sps.check()) throw NgException("Cut must be SurfacePrimitive in constructor of SplineSurface!"); auto sp = dynamic_cast(sps()->GetSolid()->GetPrimitive()); @@ -294,21 +293,21 @@ DLL_HEADER void ExportCSG(py::module &m) if(!primitive) throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!"); return make_shared(primitive,acuts); - }),bp::default_call_policies(),(bp::arg("base"), bp::arg("cuts")=bp::list()))) + }),(py::arg("base"), py::arg("cuts")=py::list())) .def("AddPoint", FunctionPointer ([] (SplineSurface & self, double x, double y, double z, bool hpref) { self.AppendPoint(Point<3>(x,y,z),hpref); return self.GetNP()-1; }), - (bp::arg("self"),bp::arg("x"),bp::arg("y"),bp::arg("z"),bp::arg("hpref")=false)) + py::arg("x"),py::arg("y"),py::arg("z"),py::arg("hpref")=false) .def("AddSegment", FunctionPointer ([] (SplineSurface & self, int i1, int i2, string bcname, double maxh) { auto str = new string(bcname); self.AppendSegment(new LineSeg<3>(self.GetPoint(i1),self.GetPoint(i2)),str,maxh); }), - (bp::arg("self"),bp::arg("pnt1"),bp::arg("pnt2"),bp::arg("bcname")="default", bp::arg("maxh")=-1.)) + py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.) ; py::class_> (m, "Solid") @@ -496,7 +495,7 @@ DLL_HEADER void ExportCSG(py::module &m) for(auto p : surf->GetPoints()) self.AddUserPoint(p); }), - (bp::arg("self"), bp::arg("SplineSurface"))) + py::arg("SplineSurface")) .def("CloseSurfaces", FunctionPointer ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, py::list aslices ) From 8a1cf75c5e7f7e30c68768f2aead51fc360857aa Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 9 Nov 2016 16:10:04 +0100 Subject: [PATCH 9/9] fix init of splinesurface --- external_dependencies/pybind11 | 2 +- libsrc/csg/python_csg.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/external_dependencies/pybind11 b/external_dependencies/pybind11 index 2b92a491..030d10e8 160000 --- a/external_dependencies/pybind11 +++ b/external_dependencies/pybind11 @@ -1 +1 @@ -Subproject commit 2b92a49115657712c7dd924248df2286e6143b8d +Subproject commit 030d10e826b87f8cdf0816aa36b9a515fb7d064d diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 9cd9390d..afc4836a 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -275,7 +275,7 @@ DLL_HEADER void ExportCSG(py::module &m) py::class_> (m, "SplineSurface", "A surface for co dim 2 integrals on the splines") - .def("__init__", FunctionPointer ([](const shared_ptr base, py::list cuts) + .def("__init__", FunctionPointer ([](SplineSurface* instance, shared_ptr base, py::list cuts) { auto primitive = dynamic_cast (base->GetSolid()->GetPrimitive()); auto acuts = new Array(); @@ -292,8 +292,8 @@ DLL_HEADER void ExportCSG(py::module &m) } if(!primitive) throw NgException("Base is not a SurfacePrimitive in constructor of SplineSurface!"); - return make_shared(primitive,acuts); - }),(py::arg("base"), py::arg("cuts")=py::list())) + new (instance) SplineSurface(primitive,acuts); + }),py::arg("base"), py::arg("cuts")=py::list()) .def("AddPoint", FunctionPointer ([] (SplineSurface & self, double x, double y, double z, bool hpref) {