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/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/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/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/edgeflw.cpp b/libsrc/csg/edgeflw.cpp index 5a5e22a6..6afa4957 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,26 @@ namespace netgen layer, mesh); } - + 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,*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); + } + + } + + } /* // not available ... @@ -1157,7 +1176,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; diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index 4b8f1713..ad67151f 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"); @@ -338,7 +344,6 @@ namespace netgen FaceDescriptor & fd = mesh.GetFaceDescriptor(k); fd.SetBCName ( mesh.GetBCNamePtr ( fd.BCProperty() - 1 ) ); } - //!! diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 60db9ee4..afc4836a 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -254,6 +254,61 @@ DLL_HEADER void ExportCSG(py::module &m) })) ; + py::class_,shared_ptr>> (m,"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])); + })) + ; + + py::class_> (m, "SplineSurface", + "A surface for co dim 2 integrals on the splines") + .def("__init__", FunctionPointer ([](SplineSurface* instance, 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]); + 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!"); + 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) + { + self.AppendPoint(Point<3>(x,y,z),hpref); + return self.GetNP()-1; + }), + 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); + }), + py::arg("pnt1"),py::arg("pnt2"),py::arg("bcname")="default", py::arg("maxh")=-1.) + ; py::class_> (m, "Solid") .def ("__str__", &ToString) @@ -427,9 +482,21 @@ DLL_HEADER void ExportCSG(py::module &m) }), py::arg("surface"), py::arg("solid") ) - + .def("AddSplineSurface", FunctionPointer + ([] (CSGeometry & self, shared_ptr surf) + { + auto cuttings = surf->CreateCuttingSurfaces(); + auto spsol = make_shared(new Solid(&*surf)); + 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); + }), + py::arg("SplineSurface")) - .def("CloseSurfaces", FunctionPointer ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, py::list aslices ) { diff --git a/libsrc/csg/splinesurface.cpp b/libsrc/csg/splinesurface.cpp new file mode 100644 index 00000000..d66c8bff --- /dev/null +++ b/libsrc/csg/splinesurface.cpp @@ -0,0 +1,67 @@ + +#include + +namespace netgen +{ +void SplineSurface :: AppendPoint(const Point<3> & p, const double reffac, const bool hpref) +{ + auto pp = Point<3>(p); + geompoints.Append(GeomPoint<3>(pp,reffac)); + geompoints.Last().hpref = hpref; +} + + void SplineSurface :: AppendSegment(SplineSeg<3>* spline, string* bcname, double amaxh) + { + splines.Append(spline); + bcnames.Append(bcname); + maxh.Append(amaxh); + } + + string* SplineSurface :: GetBCNameOf (Point<3> p1, Point<3> p2) const + { + + double eps = 1e-5; + for(int i=0; i(splines[i]->GetPoint(0)); + Project(pp1); + auto pp2 = Point<3>(splines[i]->GetPoint(1)); + Project(pp2); + if (((pp1-p1).Length()* SplineSurface :: CreateCuttingSurfaces() const + { + auto cuttings = new Array(); + for (auto cut : *cuts) + cuttings->Append(cut); + for(int i = 0; i*>(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 cuttings; + } + + void SplineSurface :: Print(ostream & str) const +{ + str << "SplineSurface " << endl; +} + +} 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..1f68ae73 --- /dev/null +++ b/libsrc/csg/splinesurface.hpp @@ -0,0 +1,79 @@ +#ifndef FILE_SPLINESURFACE +#define FILE_SPLINESURFACE + + +namespace netgen +{ + class SplineSurface : public OneSurfacePrimitive + { + protected: + Array> geompoints; + Array*> splines; + Array bcnames; + Array maxh; + OneSurfacePrimitive* baseprimitive; + Array* cuts; + + public: + 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]; } + 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, 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; + 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/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index 7ffb837d..ca272439 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)].edgenr; + 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; if (mesh->GetDimension() == 2) ret.mat = mesh->GetBCNamePtr(el.si-1); else 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..d7548127 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, @@ -519,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, @@ -565,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, @@ -584,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, @@ -682,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, @@ -691,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, diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index f190e683..2423c03e 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.edgenr-1]); + else + seg.SetBCName(0); + } + } + } + if (strcmp (str, "singular_points") == 0) { infile >> n; @@ -1359,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); @@ -5714,6 +5764,48 @@ 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()) + return defaultstring; + + 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..d09b7c98 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,12 @@ 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; + size_t GetNCD2Names() const { return cd2names.Size(); } + string * GetBCNamePtr (int bcnr) const { return bcnr < bcnames.Size() ? bcnames[bcnr] : nullptr; } 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 2dadf981..b4476e00 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -2419,10 +2419,6 @@ namespace netgen } - - - - Identifications :: Identifications (Mesh & amesh) : mesh(amesh) { diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 522eb8c8..fb6e66f0 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; } @@ -826,7 +828,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 diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 7c27092a..d1d3fab4 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -229,6 +229,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) li.append (py::cast(self.surfnr2)); return li; })) + .def_property_readonly("index", FunctionPointer([](const Segment &self) -> size_t + { + return self.edgenr; + })) ; @@ -547,6 +551,12 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) GenerateBoundaryLayer (self, blp); } )) + + .def ("Scale", FunctionPointer([](Mesh & self, double factor) + { + for(auto i = 0; i