From dc15e50956ab802ee98f87438d2bca99a636d764 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Sun, 31 May 2020 21:58:21 +0200 Subject: [PATCH 01/17] Added glueing to OCC interface, geom.Glue() from Python --- libsrc/occ/occgeom.cpp | 81 +++++++++++++++++++++++++++++++++++++++ libsrc/occ/occgeom.hpp | 3 +- libsrc/occ/python_occ.cpp | 1 + 3 files changed, 84 insertions(+), 1 deletion(-) diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 037022f9..5328bd03 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -169,8 +169,89 @@ namespace netgen } + void OCCGeometry :: GlueGeometry() + { + PrintMessage(1, "OCC Glue Geometry"); + /* + // + BRep_Builder builder; + TopoDS_Shape my_fuse; + int cnt = 0; + for (TopExp_Explorer exp_solid(shape, TopAbs_SOLID); exp_solid.More(); exp_solid.Next()) + { + cout << "cnt = " << cnt << endl; + if (cnt == 0) + my_fuse = exp_solid.Current(); + else + // my_fuse = BRepAlgoAPI_Fuse (my_fuse, exp_solid.Current()); + my_fuse = QANewModTopOpe_Glue::QANewModTopOpe_Glue(my_fuse, exp_solid.Current()); + cnt++; + } + cout << "remove" << endl; + // for (int i = 1; i <= somap.Size(); i++) + // builder.Remove (shape, somap(i)); + cout << "now add" << endl; + // builder.Add (shape, my_fuse); + shape = my_fuse; + cout << "build fmap" << endl; + BuildFMap(); + */ + // from + // https://www.opencascade.com/doc/occt-7.4.0/overview/html/occt_user_guides__boolean_operations.html + BOPAlgo_Builder aBuilder; + + // Setting arguments + TopTools_ListOfShape aLSObjects; + for (TopExp_Explorer exp_solid(shape, TopAbs_SOLID); exp_solid.More(); exp_solid.Next()) + aLSObjects.Append (exp_solid.Current()); + aBuilder.SetArguments(aLSObjects); + + // Setting options for GF + // Set parallel processing mode (default is false) + // Standard_Boolean bRunParallel = Standard_True; + // aBuilder.SetRunParallel(bRunParallel); + + // Set Fuzzy value (default is Precision::Confusion()) + // Standard_Real aFuzzyValue = 1.e-5; + // aBuilder.SetFuzzyValue(aFuzzyValue); + + // Set safe processing mode (default is false) + // Standard_Boolean bSafeMode = Standard_True; + // aBuilder.SetNonDestructive(bSafeMode); + + // Set Gluing mode for coinciding arguments (default is off) + // BOPAlgo_GlueEnum aGlue = BOPAlgo_GlueShift; + // aBuilder.SetGlue(aGlue); + + // Disabling/Enabling the check for inverted solids (default is true) + // Standard Boolean bCheckInverted = Standard_False; + // aBuilder.SetCheckInverted(bCheckInverted); + + // Set OBB usage (default is false) + // Standard_Boolean bUseOBB = Standard_True; + // aBuilder.SetUseOBB(buseobb); + + // Perform the operation + aBuilder.Perform(); + // Check for the errors + if (aBuilder.HasErrors()) + { + cout << "builder has errors" << endl; + return; + } + // Check for the warnings + if (aBuilder.HasWarnings()) + { + // treatment of the warnings + ; + } + // result of the operation + shape = aBuilder.Shape(); + BuildFMap(); + } + void OCCGeometry :: HealGeometry () { int nrc = 0, nrcs = 0, diff --git a/libsrc/occ/occgeom.hpp b/libsrc/occ/occgeom.hpp index 372979ad..447d212f 100644 --- a/libsrc/occ/occgeom.hpp +++ b/libsrc/occ/occgeom.hpp @@ -78,7 +78,7 @@ #include "Bnd_Box.hxx" #include "ShapeAnalysis.hxx" #include "ShapeBuild_ReShape.hxx" - +#include "BOPAlgo_Builder.hxx" // Philippose - 29/01/2009 // OpenCascade XDE Support @@ -343,6 +343,7 @@ namespace netgen void MakeSolid(); void HealGeometry(); + void GlueGeometry(); // Philippose - 15/01/2009 // Sets the maximum mesh size for a given face diff --git a/libsrc/occ/python_occ.cpp b/libsrc/occ/python_occ.cpp index 837ecdb0..a7508c31 100644 --- a/libsrc/occ/python_occ.cpp +++ b/libsrc/occ/python_occ.cpp @@ -66,6 +66,7 @@ DLL_HEADER void ExportNgOCC(py::module &m) }), py::arg("filename"), "Load OCC geometry from step, brep or iges file") .def(NGSPickle()) + .def("Glue", &OCCGeometry::GlueGeometry) .def("Heal",[](OCCGeometry & self, double tolerance, bool fixsmalledges, bool fixspotstripfaces, bool sewfaces, bool makesolids, bool splitpartitions) { self.tolerance = tolerance; From 1d97367e30546dbdc5f5504b0be4c917681f1108 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Tue, 2 Jun 2020 08:51:51 +0200 Subject: [PATCH 02/17] check OCC-Version of HasErrors --- libsrc/occ/occgeom.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index 5328bd03..abee8107 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -236,6 +236,7 @@ namespace netgen // Perform the operation aBuilder.Perform(); // Check for the errors +#if OCC_VERSION_HEX >= 0x070000 if (aBuilder.HasErrors()) { cout << "builder has errors" << endl; @@ -247,6 +248,7 @@ namespace netgen // treatment of the warnings ; } +#endif // result of the operation shape = aBuilder.Shape(); BuildFMap(); From 9b28a2df026753a2ba6dfcdcad3a641799c211a8 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Wed, 3 Jun 2020 11:50:33 +0200 Subject: [PATCH 03/17] OCC - HasErrors() available from v7.2 --- libsrc/occ/occgeom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/occ/occgeom.cpp b/libsrc/occ/occgeom.cpp index abee8107..91dfb2cf 100644 --- a/libsrc/occ/occgeom.cpp +++ b/libsrc/occ/occgeom.cpp @@ -236,7 +236,7 @@ namespace netgen // Perform the operation aBuilder.Perform(); // Check for the errors -#if OCC_VERSION_HEX >= 0x070000 +#if OCC_VERSION_HEX >= 0x070200 if (aBuilder.HasErrors()) { cout << "builder has errors" << endl; From 94d489e183f569f3ca97c89ebe525be31c7c076f Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Wed, 3 Jun 2020 11:54:56 +0200 Subject: [PATCH 04/17] cmake - remove compiler definition of NETGEN_VERSION --- CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a40c89c5..f67a83b2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -191,7 +191,6 @@ check_include_files (dlfcn.h HAVE_DLFCN_H) if(HAVE_DLFCN_H) add_definitions(-DHAVE_DLFCN_H) endif() -add_definitions(-DNETGEN_VERSION="${NETGEN_VERSION}") include_directories(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) From efbd71c8d58e7715d947eae3636ab48c545c971f Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 3 Jun 2020 12:42:35 +0200 Subject: [PATCH 05/17] define cmake export compile commands after project --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f67a83b2..db13de81 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -35,8 +35,6 @@ option( USE_SUPERBUILD "use ccache" ON) set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_modules") -set(CMAKE_EXPORT_COMPILE_COMMANDS ON) - if(APPLE) set(INSTALL_DIR_DEFAULT /Applications/Netgen.app) else(APPLE) @@ -78,6 +76,8 @@ else() endif() endif() +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + include (${CMAKE_CURRENT_LIST_DIR}/cmake/generate_version_file.cmake) set(CPACK_PACKAGE_VERSION "${NETGEN_VERSION}") From 5bea3bb6122a8526dec8e89bbb6799c8b2ce80ea Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:39:31 +0200 Subject: [PATCH 06/17] Implement and export SplineGeometry2d::SetDomainTensorMeshing --- libsrc/geom2d/geometry2d.hpp | 15 +++++++++++++-- libsrc/geom2d/python_geom2d.cpp | 1 + 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index ff4c459c..f5841b19 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -133,7 +133,7 @@ namespace netgen NgArray materials; NgArray maxh; NgArray quadmeshing; - NgArray tensormeshing; + Array tensormeshing; NgArray layer; NgArray bcnames; double elto0 = 1.0; @@ -216,9 +216,20 @@ namespace netgen } bool GetDomainTensorMeshing ( int domnr ) { - if ( tensormeshing.Size() ) return tensormeshing[domnr-1]; + if ( tensormeshing.Size()>=domnr ) return tensormeshing[domnr-1]; else return false; } + void SetDomainTensorMeshing ( int domnr, bool tm ) + { + if ( tensormeshing.Size()(), meshingparameter_description.c_str()) + .def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing) ; } From d1c7a16d63ba82071957aa21cae74d999c8bb6c5 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:41:22 +0200 Subject: [PATCH 07/17] Do linear interpolation of corresponding edge points in SplineGeometry tensor mesh generation better results for curved domains --- libsrc/geom2d/genmesh2d.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 09ef0260..89eaa098 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -527,11 +527,16 @@ namespace netgen for (PointIndex pix = nextpi[c1], ix = 0; pix != c2; pix = nextpi[pix], ix++) + { + Point<3> px = (*mesh)[pix]; for (PointIndex piy = nextpi[c2], iy = 0; piy != c3; piy = nextpi[piy], iy++) { - Point<3> p = (*mesh)[pix] + ( (*mesh)[piy] - (*mesh)[c2] ); - pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (p , 1, FIXEDPOINT); + double lam = Dist((*mesh)[piy],(*mesh)[c2]) / Dist((*mesh)[c3],(*mesh)[c2]); + auto pix1 = pts[(nex+1)*ney+ix+1]; + auto pnew = px + lam*((*mesh)[pix1]-px); + pts[(nex+1)*(iy+1) + ix+1] = mesh -> AddPoint (pnew, 1, FIXEDPOINT); } + } for (int i = 0; i < ney; i++) for (int j = 0; j < nex; j++) From d08e2daa060da23965f68bfc4b083ae8dd271539 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:42:26 +0200 Subject: [PATCH 08/17] Do edge swapping for faces individually with tensor product meshes If the mesh contains quads, the edge swapping algorithm switches to generic improve, which introduces quads everywhere. This is not intended if one domain contains a tensor product mesh. Thus, call the optimizer for each face if mesh contains quads but mp.quad is not set. --- libsrc/meshing/meshfunc2d.cpp | 45 ++++++++++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 3 deletions(-) diff --git a/libsrc/meshing/meshfunc2d.cpp b/libsrc/meshing/meshfunc2d.cpp index 88955a36..84f5d399 100644 --- a/libsrc/meshing/meshfunc2d.cpp +++ b/libsrc/meshing/meshfunc2d.cpp @@ -20,6 +20,19 @@ namespace netgen } mesh.Compress(); + bool optimize_swap_separate_faces = false; + if(!mp.quad) + { + bool mixed = false; + ParallelFor( Range(mesh.GetNSE()), [&] (auto i) NETGEN_LAMBDA_INLINE + { + if (mesh[SurfaceElementIndex(i)].GetNP() != 3) + mixed = true; + }); + if(mixed) + optimize_swap_separate_faces = true; + } + const char * optstr = mp.optimize2d.c_str(); int optsteps = mp.optsteps2d; @@ -31,16 +44,42 @@ namespace netgen { case 's': { // topological swap - MeshOptimize2d meshopt(mesh); + MeshOptimize2d meshopt(mesh); meshopt.SetMetricWeight (mp.elsizeweight); - meshopt.EdgeSwapping (0); + + if(optimize_swap_separate_faces) + { + for(auto i : Range(1, mesh.GetNFD()+1)) + { + meshopt.SetFaceIndex(i); + meshopt.EdgeSwapping (0); + } + } + else + { + meshopt.SetFaceIndex(0); + meshopt.EdgeSwapping (0); + } break; } case 'S': { // metric swap MeshOptimize2d meshopt(mesh); meshopt.SetMetricWeight (mp.elsizeweight); - meshopt.EdgeSwapping (1); + + if(optimize_swap_separate_faces) + { + for(auto i : Range(1, mesh.GetNFD()+1)) + { + meshopt.SetFaceIndex(i); + meshopt.EdgeSwapping (1); + } + } + else + { + meshopt.SetFaceIndex(0); + meshopt.EdgeSwapping (1); + } break; } case 'm': From c0f50820cb72d6df5d6411490884dee3a800c834 Mon Sep 17 00:00:00 2001 From: Matthias Hochsteger Date: Mon, 8 Jun 2020 10:54:02 +0200 Subject: [PATCH 09/17] Add test for SetDomainTensorMeshing --- .../test_splinegeo_tensordomainmeshing.py | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 tests/pytest/test_splinegeo_tensordomainmeshing.py diff --git a/tests/pytest/test_splinegeo_tensordomainmeshing.py b/tests/pytest/test_splinegeo_tensordomainmeshing.py new file mode 100644 index 00000000..0c00ae7d --- /dev/null +++ b/tests/pytest/test_splinegeo_tensordomainmeshing.py @@ -0,0 +1,22 @@ +from netgen.geom2d import * + +def test_tensordomainmeshing(): + geo = SplineGeometry() + w = 10 + h = 0.01 + + p = [ (0, 0), (w, 0), (w, h), (0, h) ] + p = [geo.AppendPoint(*px) for px in p] + + l0 = geo.Append ( ["line", p[0], p[1]], leftdomain=1, rightdomain=0 ) + l1 = geo.Append ( ["line", p[1], p[2]], leftdomain=1, rightdomain=0) + geo.Append ( ["line", p[3], p[2]], leftdomain=0, rightdomain=1, copy=l0 ) + geo.Append ( ["line", p[0], p[3]], leftdomain=0, rightdomain=1, copy=l1 ) + + geo._SetDomainTensorMeshing(1, True) + + mesh = geo.GenerateMesh(maxh=1) + + for el in mesh.Elements2D(): + print(el.vertices) + assert len(el.vertices) == 4 From 6a834f13ac16e0b75f9a392650de0dcc90560bff Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 16 Jun 2020 13:53:36 +0200 Subject: [PATCH 10/17] fix boundary names for boundarylayer --- libsrc/meshing/boundarylayer.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index ae7c7a30..55960395 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -163,6 +163,8 @@ namespace netgen domin, domout, -1); fd.SetBCProperty(max_surface_index); mesh.AddFaceDescriptor(fd); + mesh.SetBCName(max_surface_index-1, + "mapped_" + old_fd.GetBCName()); return max_surface_index; } return last_layer_surface_index_map[si]; From ac45a5f736105d1dd0760d5a6c06d4e85ae52f45 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 16 Jun 2020 13:54:13 +0200 Subject: [PATCH 11/17] add more information to illegal bc number exception --- 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 b38aeb88..76a22ffc 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -6541,7 +6541,7 @@ namespace netgen return defaultstring; if (bcnr < 0 || bcnr >= bcnames.Size()) - throw NgException ("illegal bc-number"); + throw RangeException("Illegal bc number ", bcnr, 0, bcnames.Size()); if ( bcnames[bcnr] ) return *bcnames[bcnr]; From 3b5c346e639237f303b731b2c9116aa0c034705c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joachim=20Sch=C3=B6berl?= Date: Wed, 17 Jun 2020 19:11:17 +0200 Subject: [PATCH 12/17] proper terms --- libsrc/core/array.hpp | 8 ++++---- libsrc/general/ngarray.hpp | 20 ++++++++++---------- libsrc/interface/readtetmesh.cpp | 18 +++++++++--------- libsrc/interface/writeabaqus.cpp | 20 ++++++++++---------- libsrc/interface/writefeap.cpp | 8 ++++---- libsrc/interface/writetet.cpp | 18 +++++++++--------- libsrc/meshing/bcfunctions.cpp | 2 +- libsrc/meshing/bcfunctions.hpp | 2 +- libsrc/meshing/parallelmesh.cpp | 6 +++--- 9 files changed, 51 insertions(+), 51 deletions(-) diff --git a/libsrc/core/array.hpp b/libsrc/core/array.hpp index ccafef6b..a1161001 100644 --- a/libsrc/core/array.hpp +++ b/libsrc/core/array.hpp @@ -1284,7 +1284,7 @@ namespace ngcore /// bubble sort array template - inline void BubbleSort (FlatArray data, FlatArray slave) + inline void BubbleSort (FlatArray data, FlatArray index) { for (size_t i = 0; i < data.Size(); i++) for (size_t j = i+1; j < data.Size(); j++) @@ -1294,9 +1294,9 @@ namespace ngcore data[i] = data[j]; data[j] = hv; - S hvs = slave[i]; - slave[i] = slave[j]; - slave[j] = hvs; + S hvs = index[i]; + index[i] = index[j]; + index[j] = hvs; } } diff --git a/libsrc/general/ngarray.hpp b/libsrc/general/ngarray.hpp index f3eea9d4..fa160a8b 100644 --- a/libsrc/general/ngarray.hpp +++ b/libsrc/general/ngarray.hpp @@ -730,7 +730,7 @@ namespace netgen /// bubble sort array template - inline void BubbleSort (NgFlatArray & data, NgFlatArray & slave) + inline void BubbleSort (NgFlatArray & data, NgFlatArray & index) { for (int i = 0; i < data.Size(); i++) for (int j = i+1; j < data.Size(); j++) @@ -740,16 +740,16 @@ namespace netgen data[i] = data[j]; data[j] = hv; - S hvs = slave[i]; - slave[i] = slave[j]; - slave[j] = hvs; + S hvs = index[i]; + index[i] = index[j]; + index[j] = hvs; } } template void QuickSortRec (NgFlatArray & data, - NgFlatArray & slave, + NgFlatArray & index, int left, int right) { int i = left; @@ -764,20 +764,20 @@ namespace netgen if (i <= j) { ngcore::Swap (data[i], data[j]); - ngcore::Swap (slave[i], slave[j]); + ngcore::Swap (index[i], index[j]); i++; j--; } } while (i <= j); - if (left < j) QuickSortRec (data, slave, left, j); - if (i < right) QuickSortRec (data, slave, i, right); + if (left < j) QuickSortRec (data, index, left, j); + if (i < right) QuickSortRec (data, index, i, right); } template - void QuickSort (NgFlatArray & data, NgFlatArray & slave) + void QuickSort (NgFlatArray & data, NgFlatArray & index) { if (data.Size() > 1) - QuickSortRec (data, slave, 0, data.Size()-1); + QuickSortRec (data, index, 0, data.Size()-1); } diff --git a/libsrc/interface/readtetmesh.cpp b/libsrc/interface/readtetmesh.cpp index e061c722..46607c59 100644 --- a/libsrc/interface/readtetmesh.cpp +++ b/libsrc/interface/readtetmesh.cpp @@ -154,7 +154,7 @@ namespace netgen break; case 7: - // NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), PID: + // NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), PID: { cout << "read nodes" << endl; for(int i=0; i> dummyint >> dummyint >> dummyint; break; @@ -254,7 +254,7 @@ namespace netgen break; case 18: - // MasterEdgeID, 3 SlaveEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2) + // MasterEdgeID, 3 MinionEdgeID's, 3 TranslCode (1=dS1 2=dS2 3=dS1+dS2) for(int i=0; i> dummyint; @@ -266,7 +266,7 @@ namespace netgen break; case 19: - // FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PSlave), PID + // FaceID, EdgeID0, EdgeID1, EdgeID2, FaceType (0=Reg 1=PMaster 2=PMinion), PID { //Segment seg; int segnum_ng[3]; @@ -343,7 +343,7 @@ namespace netgen break; case 21: - // MasterFaceID, SlaveFaceID, TranslCode (1=dS1 2=dS2) + // MasterFaceID, MinionFaceID, TranslCode (1=dS1 2=dS2) { Vec<3> randomvec(-1.32834,3.82399,0.5429151); int maxtransl = -1; diff --git a/libsrc/interface/writeabaqus.cpp b/libsrc/interface/writeabaqus.cpp index 916600af..209ac89a 100644 --- a/libsrc/interface/writeabaqus.cpp +++ b/libsrc/interface/writeabaqus.cpp @@ -158,25 +158,25 @@ void WriteAbaqusFormat (const Mesh & mesh, cout << "masternode = " << masternode << " = " << mesh.Point(masternode) << endl; - NgArray slaves(3); + NgArray minions(3); for (i = 1; i <= 3; i++) { mesh.GetIdentifications().GetPairs (i, pairs); for (j = 1; j <= pairs.Size(); j++) { if (pairs.Get(j).I1() == masternode) - slaves.Elem(i) = pairs.Get(j).I2(); + minions.Elem(i) = pairs.Get(j).I2(); } - cout << "slave(" << i << ") = " << slaves.Get(i) - << " = " << mesh.Point(slaves.Get(i)) << endl; + cout << "minion(" << i << ") = " << minions.Get(i) + << " = " << mesh.Point(minions.Get(i)) << endl; } outfile << "**\n" << "*NSET,NSET=CTENODS\n" - << slaves.Get(1) << ", " - << slaves.Get(2) << ", " - << slaves.Get(3) << endl; + << minions.Get(1) << ", " + << minions.Get(2) << ", " + << minions.Get(3) << endl; outfile << "**\n" @@ -190,7 +190,7 @@ void WriteAbaqusFormat (const Mesh & mesh, << "*BOUNDARY, OP=NEW\n"; for (j = 1; j <= 3; j++) { - Vec3d v(mesh.Point(masternode), mesh.Point(slaves.Get(j))); + Vec3d v(mesh.Point(masternode), mesh.Point(minions.Get(j))); double vlen = v.Length(); int dir = 0; if (fabs (v.X()) > 0.9 * vlen) dir = 2; @@ -198,7 +198,7 @@ void WriteAbaqusFormat (const Mesh & mesh, if (fabs (v.Z()) > 0.9 * vlen) dir = 1; if (!dir) cout << "ERROR: Problem with rigid body constraints" << endl; - outfile << slaves.Get(j) << ", " << dir << ",, 0.\n"; + outfile << minions.Get(j) << ", " << dir << ",, 0.\n"; } outfile << "**\n" @@ -223,7 +223,7 @@ void WriteAbaqusFormat (const Mesh & mesh, mpc << "4" << "\n"; mpc << pairs.Get(j).I2() << "," << k << ", -1.0, "; mpc << pairs.Get(j).I1() << "," << k << ", 1.0, "; - mpc << slaves.Get(i) << "," << k << ", 1.0, "; + mpc << minions.Get(i) << "," << k << ", 1.0, "; mpc << masternode << "," << k << ", -1.0 \n"; } } diff --git a/libsrc/interface/writefeap.cpp b/libsrc/interface/writefeap.cpp index dc2574a9..84f71eef 100644 --- a/libsrc/interface/writefeap.cpp +++ b/libsrc/interface/writefeap.cpp @@ -144,11 +144,11 @@ void WriteFEAPFormat (const Mesh & mesh, // BEGIN CONTACT OUTPUT /* - int masterindex, slaveindex; + int masterindex, minionindex; cout << "Master Surface index = "; cin >> masterindex; - cout << "Slave Surface index = "; - cin >> slaveindex; + cout << "Minion Surface index = "; + cin >> minionindex; // CONTACT SURFACE 1 @@ -196,7 +196,7 @@ void WriteFEAPFormat (const Mesh & mesh, Element2d sel = mesh.SurfaceElement(i); if (invertsurf) sel.Invert(); - if (mesh.GetFaceDescriptor(sel.GetIndex ()).BCProperty() == slaveindex) + if (mesh.GetFaceDescriptor(sel.GetIndex ()).BCProperty() == minionindex) { zz++; outfile.width(14); diff --git a/libsrc/interface/writetet.cpp b/libsrc/interface/writetet.cpp index 3bc0ba5a..2dbfa440 100644 --- a/libsrc/interface/writetet.cpp +++ b/libsrc/interface/writetet.cpp @@ -427,7 +427,7 @@ namespace netgen << numedges << " " << numnodes << endl << endl; - outfile << "// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PSlave 3=CPMaster 4=CPSlave), "<< uidpid <<":\n" \ + outfile << "// NodeID, X, Y, Z, Type (0=Reg 1=PMaster 2=PMinion 3=CPMaster 4=CPMinion), "<< uidpid <<":\n" \ << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"; @@ -515,7 +515,7 @@ namespace netgen << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n" \ << n2 << "\n" \ << "\n" \ - << "// MasterNodeID, SlaveNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n" \ + << "// MasterNodeID, MinionNodeID, TranslCode (1=dS1 2=dS2 3=dS1+dS2):\n" \ << "// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^\n"; for(int i=0; i & volume_weights , NgArray & surface_weights, NgArray & segment_weights) { From d2cb67f681d8ba4635805445e89fbb5909f6ba9a Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Fri, 19 Jun 2020 17:36:48 +0200 Subject: [PATCH 13/17] fix cmake warning --- cmake/cmake_modules/FindOpenCasCade.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/cmake_modules/FindOpenCasCade.cmake b/cmake/cmake_modules/FindOpenCasCade.cmake index e82f8dac..870e8ed5 100644 --- a/cmake/cmake_modules/FindOpenCasCade.cmake +++ b/cmake/cmake_modules/FindOpenCasCade.cmake @@ -27,7 +27,7 @@ endif(WIN32) if(OCC_LIBRARY AND NOT OCC_LIBRARY_DIR) get_filename_component(OCC_LIBRARY_DIR ${OCC_LIBRARY} PATH) -endif(OCC_LIBRARY) +endif(OCC_LIBRARY AND NOT OCC_LIBRARY_DIR) if(OCC_INCLUDE_DIR) file(STRINGS ${OCC_INCLUDE_DIR}/Standard_Version.hxx OCC_MAJOR From c3441344fbf86f7c7cc3a9a7b892230b69f6a54e Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Tue, 23 Jun 2020 18:52:29 +0200 Subject: [PATCH 14/17] set material in tensorproduct mesh in 2d as well --- libsrc/geom2d/genmesh2d.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 89eaa098..c1eadf1a 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -550,6 +550,10 @@ namespace netgen mesh -> AddSurfaceElement (el); } + char* material; + geometry.GetMaterial(domnr, material); + if(material) + mesh->SetMaterial(domnr, material); } From 177ecc74594456fa66a438099c3ed3a722a3f3b4 Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Wed, 24 Jun 2020 06:41:06 +0000 Subject: [PATCH 15/17] Allow curving of mesh if boundarylayer is flat. If surfnr is larger than nr of surfaces then do linear interpolation for PointInBetween and so on. Some fixes in boundarylayer so that surface numbers are correct. --- libsrc/csg/csgeom.cpp | 1 + libsrc/meshing/boundarylayer.cpp | 132 +++++++++++++------------------ libsrc/meshing/python_mesh.cpp | 18 ++++- 3 files changed, 74 insertions(+), 77 deletions(-) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 4874f0a9..cc3e4464 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -130,6 +130,7 @@ namespace netgen Point<3> & newp, EdgePointGeomInfo & newgi) const { Point<3> hnewp = p1+secpoint*(p2-p1); + //(*testout) << "hnewp " << hnewp << " s1 " << surfi1 << " s2 " << surfi2 << endl; if (surfi1 != -1 && surfi2 != -1 && surfi1 != surfi2) { diff --git a/libsrc/meshing/boundarylayer.cpp b/libsrc/meshing/boundarylayer.cpp index 55960395..3ce5891a 100644 --- a/libsrc/meshing/boundarylayer.cpp +++ b/libsrc/meshing/boundarylayer.cpp @@ -159,7 +159,8 @@ namespace netgen auto& old_fd = mesh.GetFaceDescriptor(si); int domout = blp.outside ? old_fd.DomainOut() : blp.new_matnrs[layer-1]; int domin = blp.outside ? blp.new_matnrs[layer-1] : old_fd.DomainIn(); - FaceDescriptor fd(max_surface_index-1, + // -1 surf nr is so that curving does not do anything + FaceDescriptor fd(-1, domin, domout, -1); fd.SetBCProperty(max_surface_index); mesh.AddFaceDescriptor(fd); @@ -271,55 +272,53 @@ namespace netgen if(blp.grow_edges) for(SegmentIndex sei = 0; sei < nseg; sei++) { - PointIndex seg_p1 = mesh[sei][0]; - PointIndex seg_p2 = mesh[sei][1]; + auto& segi = mesh[sei]; // Only go in if the segment is still active, and if both its // surface index is part of the "hit-list" if(segsel.Test(sei)) { - if(blp.surfid.Contains(mesh[sei].si)) - { - // clear the bit to indicate that this segment has been processed - segsel.Clear(sei); - - // Find matching segment pair on other surface - for(SegmentIndex sej = 0; sej < nseg; sej++) + if(blp.surfid.Contains(segi.si)) { - PointIndex segpair_p1 = mesh[sej][1]; - PointIndex segpair_p2 = mesh[sej][0]; + // clear the bit to indicate that this segment has been processed + segsel.Clear(sei); - // Find the segment pair on the neighbouring surface element - // Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0] - if(segsel.Test(sej) && ((segpair_p1 == seg_p1) && (segpair_p2 == seg_p2))) + // Find matching segment pair on other surface + for(SegmentIndex sej = 0; sej < nseg; sej++) { - // clear bit to indicate that processing of this segment is done - segsel.Clear(sej); + auto& segj = mesh[sej]; + // Find the segment pair on the neighbouring surface element + // Identified by: seg1[0] = seg_pair[1] and seg1[1] = seg_pair[0] + if(segsel.Test(sej) && ((segi[0] == segj[1]) && (segi[1] == segj[0]))) + { + // clear bit to indicate that processing of this segment is done + segsel.Clear(sej); - // Only worry about those surfaces which are not in the - // boundary layer list - if(!blp.surfid.Contains(mesh[sej].si)) + // if segj is not in surfel list we nned to add quads + if(!blp.surfid.Contains(segj.si)) { SurfaceElementIndex pnt_commelem; SetInvalid(pnt_commelem); - auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segpair_p1); - auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segpair_p2); + auto pnt1_elems = meshtopo.GetVertexSurfaceElements(segj[0]); + auto pnt2_elems = meshtopo.GetVertexSurfaceElements(segj[1]); for(auto pnt1_sei : pnt1_elems) - if(mesh[pnt1_sei].GetIndex() == mesh[sej].si) + if(mesh[pnt1_sei].GetIndex() == segj.si) for(auto pnt2_sei : pnt2_elems) if(pnt1_sei == pnt2_sei) pnt_commelem = pnt1_sei; if(IsInvalid(pnt_commelem)) - throw Exception("Couldn't find element on other side for " + ToString(segpair_p1) + " to " + ToString(segpair_p2)); + throw Exception("Couldn't find element on other side for " + ToString(segj[0]) + " to " + ToString(segj[1])); const auto& commsel = mesh[pnt_commelem]; auto surfelem_vect = GetSurfaceNormal(mesh, commsel); if(blp.outside) surfelem_vect *= -1; Element2d sel(QUAD); + auto seg_p1 = segi[0]; + auto seg_p2 = segi[1]; if(blp.outside) Swap(seg_p1, seg_p2); sel[0] = seg_p1; @@ -332,7 +331,7 @@ namespace netgen { domains_to_surf_index[domains] = ++max_surface_index; domains_to_surf_index[make_tuple(max_surface_index, get<1>(domains), get<2>(domains))] = max_surface_index; - FaceDescriptor fd(max_surface_index-1, + FaceDescriptor fd(-1, get<1>(domains), get<2>(domains), -1); @@ -367,46 +366,37 @@ namespace netgen seg_2.edgenr = pi_to_edgenr[points]; seg_2.si = new_index; mesh.AddSegment(seg_2); - mesh[sej].si = new_index; } // in last layer insert new segments if(layer == blp.heights.Size()) { - Segment s1 = mesh[sei]; - Segment s2 = mesh[sej]; + Segment s1 = segi; + Segment s2 = segj; s1.edgenr = ++max_edge_nr; s2.edgenr = max_edge_nr; - bool create_it = true; - if(blp.surfid.Contains(mesh[sej].si)) - { - if(last_layer_surface_index_map.find(s1.si) != last_layer_surface_index_map.end() && - last_layer_surface_index_map.find(s2.si) != last_layer_surface_index_map.end()) - // edge already mapped - create_it = false; - s2.si = map_surface_index(s2.si); - } + if(blp.surfid.Contains(segj.si)) + s2.si = map_surface_index(segj.si); else { - s2.si = domains_to_surf_index[make_tuple(s2.si, - blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())]; + auto side_surf = domains_to_surf_index[make_tuple(s2.si, blp.new_matnrs[layer-1], mesh.GetFaceDescriptor(s2.si).DomainOut())]; + if(blp.outside) + s2.si = side_surf; + else + segj.si = side_surf; } s1.si = map_surface_index(s1.si); - if(create_it) - { - mesh.AddSegment(s1); - mesh.AddSegment(s2); - } + s1.surfnr1 = s1.surfnr2 = s2.surfnr1 = s2.surfnr2 = -1; + mesh.AddSegment(s1); + mesh.AddSegment(s2); + } + // segi[0] = mapto[segi[0]] not working somehow? + mesh[sei][0] = mapto[segi[0]]; + mesh[sei][1] = mapto[segi[1]]; + mesh[sej][0] = mapto[segj[0]]; + mesh[sej][1] = mapto[segj[1]]; } - - // remap the segments to the new points - mesh[sei][0] = mapto[mesh[sei][0]]; - mesh[sei][1] = mapto[mesh[sei][1]]; - mesh[sej][1] = mapto[mesh[sej][1]]; - mesh[sej][0] = mapto[mesh[sej][0]]; - } - } } else { @@ -457,9 +447,9 @@ namespace netgen { for(SurfaceElementIndex si = 0; si < nse; si++) { - if(blp.surfid.Contains(mesh[si].GetIndex())) + const auto& sel = mesh[si]; + if(blp.surfid.Contains(sel.GetIndex())) { - const auto& sel = mesh[si]; Element2d newel = sel; newel.SetIndex(map_surface_index(sel.GetIndex())); mesh.AddSurfaceElement(newel); @@ -548,36 +538,28 @@ namespace netgen for(SurfaceElementIndex sei : Range(nse)) { auto& sel = mesh[sei]; - bool to_move = blp.surfid.Contains(sel.GetIndex()); - if(blp.domains.Size()) + if(!blp.surfid.Contains(sel.GetIndex())) { - if(blp.outside) - to_move |= blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()]; - else - to_move |= !blp.domains[mesh.GetFaceDescriptor(sel.GetIndex()).DomainIn()]; - } - - if(to_move) - { - for(auto& pnum : sel.PNums()) - // Check (Doublecheck) if the corresponding point has a - // copy available for remapping - if(mapto[pnum].IsValid()) - // Map the surface elements to the new points - pnum = mapto[pnum]; + const auto& fd = mesh.GetFaceDescriptor(sel.GetIndex()); + if(blp.outside && + (!blp.domains[fd.DomainIn()] && !blp.domains[fd.DomainOut()])) + continue; + if(!blp.outside && + (blp.domains[fd.DomainIn()] || blp.domains[fd.DomainOut()])) + continue; } + for(auto& pnum : sel.PNums()) + if(mapto[pnum].IsValid()) + pnum = mapto[pnum]; } for(ElementIndex ei : Range(ne)) { auto& el = mesh[ei]; - bool to_move = blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()]; - if(blp.domains.Size() == 0 || to_move) + // only move the elements on the correct side + if(blp.outside ? blp.domains[el.GetIndex()] : !blp.domains[el.GetIndex()]) for(auto& pnum : el.PNums()) - // Check (Doublecheck) if the corresponding point has a - // copy available for remapping if(mapto[pnum].IsValid()) - // Map the volume elements to the new points pnum = mapto[pnum]; } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index fdfb2c40..60c719e7 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -949,8 +949,22 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) { regex pattern(*get_if(&boundary)); for(int i = 1; i<=self.GetNFD(); i++) - if(regex_match(self.GetFaceDescriptor(i).GetBCName(), pattern)) - blp.surfid.Append(i); + { + auto& fd = self.GetFaceDescriptor(i); + if(regex_match(fd.GetBCName(), pattern)) + { + auto dom_pattern = get_if(&domain); + // only add if adjacent to domain + if(dom_pattern) + { + regex pattern(*dom_pattern); + if(regex_match(self.GetMaterial(fd.DomainIn()), pattern) || (fd.DomainOut() > 0 ? regex_match(self.GetMaterial(fd.DomainOut()), pattern) : false)) + blp.surfid.Append(i); + } + else + blp.surfid.Append(i); + } + } } if(double* pthickness = get_if(&thickness); pthickness) From 1a619841b2b9da67f1b851fd1e004cb0f34f986c Mon Sep 17 00:00:00 2001 From: Michael Neunteufel Date: Wed, 24 Jun 2020 06:41:55 +0000 Subject: [PATCH 16/17] Surface geom --- libsrc/meshing/CMakeLists.txt | 4 +- libsrc/meshing/curvedelems.cpp | 30 +-- libsrc/meshing/meshing.hpp | 2 + libsrc/meshing/python_mesh.cpp | 39 +++ libsrc/meshing/surfacegeom.cpp | 419 +++++++++++++++++++++++++++++++++ libsrc/meshing/surfacegeom.hpp | 70 ++++++ 6 files changed, 547 insertions(+), 17 deletions(-) create mode 100644 libsrc/meshing/surfacegeom.cpp create mode 100644 libsrc/meshing/surfacegeom.hpp diff --git a/libsrc/meshing/CMakeLists.txt b/libsrc/meshing/CMakeLists.txt index 47b99e79..9bf45a89 100644 --- a/libsrc/meshing/CMakeLists.txt +++ b/libsrc/meshing/CMakeLists.txt @@ -12,7 +12,7 @@ add_library(mesh ${NG_LIB_TYPE} smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp topology.cpp triarls.cpp validate.cpp bcfunctions.cpp parallelmesh.cpp paralleltop.cpp paralleltop.hpp basegeom.cpp - python_mesh.cpp hexarls.cpp + python_mesh.cpp hexarls.cpp surfacegeom.cpp ../../ng/onetcl.cpp ${mesh_object_libs} ) @@ -37,6 +37,6 @@ install(FILES localh.hpp meshclass.hpp meshfunc.hpp meshing2.hpp meshing3.hpp meshing.hpp meshtool.hpp meshtype.hpp msghandler.hpp paralleltop.hpp ruler2.hpp ruler3.hpp specials.hpp topology.hpp validate.hpp - python_mesh.hpp + python_mesh.hpp surfacegeom.hpp DESTINATION ${NG_INSTALL_DIR_INCLUDE}/meshing COMPONENT netgen_devel ) diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 6523c056..7ea1086b 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -748,21 +748,19 @@ namespace netgen for (int i2 = 0; i2 < edgenrs.Size(); i2++) { - // PointIndex pi1 = el[edges[i2][0]]; - // PointIndex pi2 = el[edges[i2][1]]; - - // bool swap = pi1 > pi2; - - // Point<3> p1 = mesh[pi1]; - // Point<3> p2 = mesh[pi2]; - - // int order1 = edgeorder[edgenrs[i2]]; - // int ndof = max (0, order1-1); - - surfnr[edgenrs[i2]] = mesh.GetFaceDescriptor(el.GetIndex()).SurfNr(); - gi0[edgenrs[i2]] = el.GeomInfoPi(edges[i2][0]+1); - gi1[edgenrs[i2]] = el.GeomInfoPi(edges[i2][1]+1); - } + auto enr = edgenrs[i2]; + surfnr[enr] = mesh.GetFaceDescriptor(el.GetIndex()).SurfNr(); + if (el[edges[i2][0]] < el[edges[i2][1]]) + { + gi0[enr] = el.GeomInfoPi(edges[i2][0]+1); + gi1[enr] = el.GeomInfoPi(edges[i2][1]+1); + } + else + { + gi1[enr] = el.GeomInfoPi(edges[i2][0]+1); + gi0[enr] = el.GeomInfoPi(edges[i2][1]+1); + } + } } @@ -1303,6 +1301,8 @@ namespace netgen SurfaceElementIndex sei = top.GetFace2SurfaceElement (f+1)-1; if (sei != SurfaceElementIndex(-1)) { PointGeomInfo gi = mesh[sei].GeomInfoPi(1); + gi.u = 1.0/3.0*(mesh[sei].GeomInfoPi(1).u+mesh[sei].GeomInfoPi(2).u+mesh[sei].GeomInfoPi(3).u); + gi.v = 1.0/3.0*(mesh[sei].GeomInfoPi(1).v+mesh[sei].GeomInfoPi(2).v+mesh[sei].GeomInfoPi(3).v); geo.ProjectPointGI(surfnr[facenr], pp, gi); } else diff --git a/libsrc/meshing/meshing.hpp b/libsrc/meshing/meshing.hpp index 35cb9add..253f2f93 100644 --- a/libsrc/meshing/meshing.hpp +++ b/libsrc/meshing/meshing.hpp @@ -57,10 +57,12 @@ namespace netgen #include "hprefinement.hpp" #include "boundarylayer.hpp" #include "specials.hpp" + } #include "validate.hpp" #include "basegeom.hpp" +#include "surfacegeom.hpp" #ifdef PARALLEL #include "paralleltop.hpp" diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index fdfb2c40..6584878e 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1132,6 +1132,45 @@ grow_edges : bool = False m.def("ReadCGNSFile", &ReadCGNSFile, py::arg("filename"), py::arg("base")=1, "Read mesh and solution vectors from CGNS file"); + + py::class_> (m, "SurfaceGeometry") + .def(py::init<>()) + .def(py::init([](py::object pyfunc) + { + std::function (Point<2>)> func = [pyfunc](Point<2> p) + { + py::gil_scoped_acquire aq; + py::tuple pyres = py::extract(pyfunc(p[0],p[1],0.0)) (); + return Vec<3>(py::extract(pyres[0])(),py::extract(pyres[1])(),py::extract(pyres[2])()); + }; + auto geo = make_shared(func); + return geo; + }), py::arg("mapping")) + .def(NGSPickle()) + .def("GenerateMesh", [](shared_ptr geo, + bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames) + { + if (py::len(py_bbbpts) != py::len(py_bbbnames)) + throw Exception("In SurfaceGeometry::GenerateMesh bbbpts and bbbnames do not have same lengths."); + Array> bbbpts(py::len(py_bbbpts)); + Array bbbname(py::len(py_bbbpts)); + for(int i = 0; i(py_bbbpts[i])(); + bbbpts[i] = Point<3>(py::extract(pnt[0])(),py::extract(pnt[1])(),py::extract(pnt[2])()); + bbbname[i] = py::extract(py_bbbnames[i])(); + } + auto mesh = make_shared(); + SetGlobalMesh (mesh); + mesh->SetGeometry(geo); + ng_geometry = geo; + auto result = geo->GenerateMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname); + if(result != 0) + throw Exception("SurfaceGeometry: Meshing failed!"); + return mesh; + }, py::arg("quads")=true, py::arg("nx")=10, py::arg("ny")=10, py::arg("flip_triangles")=false, py::arg("bbbpts")=py::list(), py::arg("bbbnames")=py::list()) + ; + ; } PYBIND11_MODULE(libmesh, m) { diff --git a/libsrc/meshing/surfacegeom.cpp b/libsrc/meshing/surfacegeom.cpp new file mode 100644 index 00000000..40d272a8 --- /dev/null +++ b/libsrc/meshing/surfacegeom.cpp @@ -0,0 +1,419 @@ +/* *************************************************************************/ +/* File: surfacegeom.cpp */ +/* Author: Michael Neunteufel */ +/* Date: Jun. 2020 */ +/* *************************************************************************/ + +#include + +namespace netgen +{ + SurfaceGeometry :: SurfaceGeometry() + { + //identity + func = [](Point<2> p) { return Vec<3>(p[0],p[1],0.0); }; + } + + SurfaceGeometry :: SurfaceGeometry(function(Point<2>)> _func) : func(_func) + { + ; + } + + SurfaceGeometry :: SurfaceGeometry(const SurfaceGeometry& geom) : func(geom.func), eps(geom.eps) + { + ; + } + + void SurfaceGeometry :: CalcHesse(double u, double v, Vec<3>& f_uu, Vec<3>& f_vv, Vec<3>& f_uv) const + { + Point<2> p = Point<2>(u,v); + double pr = p[0]+eps; + double pl = p[0]-eps; + double prr = p[0]+2*eps; + double pll = p[0]-2*eps; + + auto dr = GetTangentVectors( pr, v ); + auto dl = GetTangentVectors( pl, v ); + auto drr = GetTangentVectors( prr, v ); + auto dll = GetTangentVectors( pll, v ); + + f_uu = (1.0/(12.0*eps)) * (8.0*dr[0]-8.0*dl[0]-drr[0]+dll[0]); + f_uv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]); + + pr = p[1]+eps; + pl = p[1]-eps; + prr = p[1]+2*eps; + pll = p[1]-2*eps; + + dr = GetTangentVectors(u, pr); + dl = GetTangentVectors(u, pl); + drr = GetTangentVectors(u, prr); + dll = GetTangentVectors(u, pll); + + f_vv = (1.0/(12.0*eps)) * (8.0*dr[1]-8.0*dl[1]-drr[1]+dll[1]); + } + + Array> SurfaceGeometry :: GetTangentVectors(double u, double v) const + { + Array> tang(2); + + Point<2> pru = Point<2>(u+eps,v); + Point<2> plu = Point<2>(u-eps,v); + Point<2> prru = Point<2>(u+2*eps,v); + Point<2> pllu = Point<2>(u-2*eps,v); + + Point<2> prv = Point<2>(u,v+eps); + Point<2> plv = Point<2>(u,v-eps); + Point<2> prrv = Point<2>(u,v+2*eps); + Point<2> pllv = Point<2>(u,v-2*eps); + + + tang[0] = 1/(12.0*eps)*( 8.0*func(pru) - 8.0*func(plu) - func(prru) + func(pllu) ); + tang[1] = 1/(12.0*eps)*( 8.0*func(prv) - 8.0*func(plv) - func(prrv) + func(pllv) ); + + return tang; + } + + Vec<3> SurfaceGeometry :: GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const + { + Array> tang = GetTangentVectors(gi->u, gi->v); + auto normal = Cross(tang[0], tang[1]); + return Cross(tang[0], tang[1]); + } + + + PointGeomInfo SurfaceGeometry :: ProjectPoint(int surfind, Point<3> & p) const + { + throw Exception("In SurfaceGeometry::ProjectPoint"); + } + + void SurfaceGeometry :: ProjectPointEdge (int surfind, int surfind2, Point<3> & p, + EdgePointGeomInfo* gi) const + { + if (gi == nullptr) + throw Exception("In SurfaceGeometry::ProjectPointEdge: gi is nullptr"); + throw Exception("In SurfaceGeometry::ProjectPointEdge: not implemented"); + } + + bool SurfaceGeometry :: ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const + { + Array> tangs; + Vec<3> diff, f_uu, f_vv, f_uv; + Vec<2> r, dx; + double norm_r, det, energy=0.0, new_energy=0.0, alpha=2.0,u=0.0,v=0.0; + Mat<2,2> mat, inv; + int num=0, maxit=20; + double damping=0.2; + + + //Solve minimization problem + // argmin_(u,v) 0.5*\| f(u,v)-p\|^2 + //via Neton's method: + // F(u,v) = ( (f(u,v)-p)*f_u(u,v), (f(u,v)-p)*f_v(u,v))^T = (0,0)^T + //Stiffness matrix + // F'(u,v) = ( f_u*f_u + (f-p)*f_uu, f_v*f_u + (f-p)*f_uv, f_v*f_u + (f-p)*f_uv, f_v*f_v + (f-p)*f_vv ) + do + { + num++; + tangs = GetTangentVectors(gi.u, gi.v); + diff = func(Point<2>(gi.u, gi.v)) - Vec<3>(p); + energy = diff.Length2(); + r = Vec<2>( diff*tangs[0], diff*tangs[1] ); + norm_r = r.Length2(); + + CalcHesse(gi.u, gi.v, f_uu, f_vv, f_uv); + + + mat(0,0) = tangs[0]*tangs[0] + diff*f_uu; + mat(1,0) = mat(0,1) = tangs[0]*tangs[1]+diff*f_uv; + mat(1,1) = tangs[1]*tangs[1]+diff*f_vv; + + CalcInverse(mat,inv); + + dx = inv*r; + + //Linesearch + alpha = 2.0; + do + { + alpha /= 2.0; + u = gi.u - min(1.0,alpha*damping*num)*dx[0]; + v = gi.v - min(1.0,alpha*damping*num)*dx[1]; + + diff = func(Point<2>(u, v)) - Vec<3>(p); + new_energy = diff.Length2(); + } + while (alpha > 1e-10 && new_energy > energy+1e-14); + if (alpha <= 1e-10) + throw Exception("In SurfaceGeometry::ProjectPointGI: Linesearch min alpha reached!"); + gi.u = u; + gi.v = v; + + + } + while ( norm_r > 1e-12 && num < maxit); + + //Stay in reference domain [0,1]^2 + if (gi.u < 0 || gi.u > 1 || gi.v < 0 || gi.v > 1) + { + cout << "Warning: Projected point outside [0,1]^2: u=" << gi.u << ",v=" << gi.v <<". Setting back." << endl; + gi.u = min(max(gi.u,0.0),1.0); + gi.v = min(max(gi.v,0.0),1.0); + } + + p = Point<3>(func(Point<2>(gi.u,gi.v))); + + if (num == maxit) + { + //cout << "In SurfaceGeometry::ProjectPointGI: Newton did not converge" << endl; + throw Exception("In SurfaceGeometry::ProjectPointGI: Newton did not converge"); + } + return true; + } + + bool SurfaceGeometry :: CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const + { + throw Exception("In SurfaceGeometry::CalcPointGeomInfo: not implemented"); + return false; + } + + void SurfaceGeometry :: PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi1, int surfi2, const EdgePointGeomInfo & ap1, const EdgePointGeomInfo & ap2, Point<3> & newp, EdgePointGeomInfo & newgi) const + { + newp = p1+secpoint*(p2-p1); + + PointGeomInfo pgi; + pgi.u = ap1.u+secpoint*(ap2.u-ap1.u); + pgi.v = ap1.v+secpoint*(ap2.v-ap1.v); + + ProjectPointGI(surfi1, newp, pgi); + + newgi.u = pgi.u; + newgi.v = pgi.v; + newgi.edgenr = ap1.edgenr; + newgi.body = -1; + newgi.dist = -1.0; + } + + void SurfaceGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, + int surfi, + const PointGeomInfo & gi1, + const PointGeomInfo & gi2, + Point<3> & newp, PointGeomInfo & newgi) const + { + newp = p1+secpoint*(p2-p1); + + newgi.u = gi1.u+secpoint*(gi2.u-gi1.u); + newgi.v = gi1.v+secpoint*(gi2.v-gi1.v); + newgi.trignum = -1; + + ProjectPointGI(surfi, newp, newgi); + } + + int SurfaceGeometry :: GenerateMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames) + { + mesh->SetDimension(3); + + Array found(bbbpts.Size()); + found = false; + Array indbbbpts(bbbpts.Size()); + + + Array pids; + Array pgis; + for(int i=0; i <= ny; i++) + for(int j=0; j <= nx; j++) + { + PointGeomInfo pgi; + pgi.trignum = -1; + pgi.u = double(j)/nx; + pgi.v = double(i)/ny; + + Point<3> pnt = Point<3>(func(Point<2>(pgi.u,pgi.v))); + pids.Append(mesh->AddPoint(pnt)); + pgis.Append(pgi); + + for (int k = 0; k < bbbpts.Size(); k++) + { + auto diff = pnt - bbbpts[k]; + if(diff.Length2() < 1e-14) + { + found[k] = true; + indbbbpts[k] = pids[pids.Size()-1]; + } + } + } + + for (bool f : found) + if (!f) + throw Exception("In SurfaceGeometry :: GenerateMesh: bbbpts not resolved in mesh."); + + FaceDescriptor fd; + fd.SetSurfNr(1); + fd.SetDomainIn(1); + fd.SetDomainOut(0); + fd.SetBCProperty(1); + mesh->AddFaceDescriptor(fd); + + + for(int i=0; i < ny; i++) + { + for(int j=0; j < nx; j++) + { + int base = i * (nx+1) + j; + if (quads) + { + int pnum[4] = {base,base+1,base+nx+2,base+nx+1}; + Element2d el = Element2d(QUAD); + for (int i = 0; i < 4; i++) + { + el[i] = pids[pnum[i]]; + el.GeomInfoPi(i+1) = pgis[pnum[i]]; + } + el.SetIndex(1); + + mesh->AddSurfaceElement(el); + } + else + { + Array pnum1(3); + Array pnum2(3); + if (flip_triangles) + { + pnum1[0] = base; + pnum1[1] = base+1; + pnum1[2] = base+nx+2; + pnum2[0] = base; + pnum2[1] = base+nx+2; + pnum2[2] = base+nx+1; + } + else + { + pnum1[0] = base; + pnum1[1] = base+1; + pnum1[2] = base+nx+1; + pnum2[0] = base+1; + pnum2[1] = base+nx+2; + pnum2[2] = base+nx+1; + } + + Element2d el = Element2d(TRIG); + for (int i = 0; i < 3; i++) + { + el[i] = pids[pnum1[i]]; + el.GeomInfoPi(i+1) = pgis[pnum1[i]]; + } + el.SetIndex(1); + + mesh->AddSurfaceElement(el); + for (int i = 0; i < 3; i++) + { + el[i] = pids[pnum2[i]]; + el.GeomInfoPi(i+1) = pgis[pnum2[i]]; + } + mesh->AddSurfaceElement(el); + } + } + } + + Segment seg; + seg.si = 1; + seg.edgenr = 1; + seg.epgeominfo[0].edgenr = 1; + seg.epgeominfo[1].edgenr = 1; + // needed for codim2 in 3d + seg.edgenr = 1; + for(int i=0; i < nx; i++) + { + seg[0] = pids[i]; + seg[1] = pids[i+1]; + + seg.geominfo[0] = pgis[i]; + seg.geominfo[1] = pgis[i+1]; + seg.epgeominfo[0].u = pgis[i].u; + seg.epgeominfo[0].v = pgis[i].v; + seg.epgeominfo[0].edgenr = seg.edgenr; + seg.epgeominfo[1].u = pgis[i+1].u; + seg.epgeominfo[1].v = pgis[i+1].v; + seg.epgeominfo[1].edgenr = seg.edgenr; + + mesh->AddSegment(seg); + } + + seg.si = 2; + seg.edgenr = 2; + for(int i=0; iAddSegment(seg); + } + + seg.si = 3; + seg.edgenr = 3; + for(int i=0; iAddSegment(seg); + } + + seg.si = 4; + seg.edgenr = 4; + for(int i=0; iAddSegment(seg); + } + + mesh->SetCD2Name(1, "bottom"); + mesh->SetCD2Name(2, "right"); + mesh->SetCD2Name(3, "top"); + mesh->SetCD2Name(4, "left"); + + for (int i = 0; i < bbbpts.Size(); i++) + { + Element0d el; + el.pnum = indbbbpts[i]; + el.index = i+1; + mesh->pointelements.Append(el); + mesh->SetCD3Name(i+1, bbbnames[i]); + } + + mesh->Compress(); + mesh->UpdateTopology(); + + return 0; + } + +}; diff --git a/libsrc/meshing/surfacegeom.hpp b/libsrc/meshing/surfacegeom.hpp new file mode 100644 index 00000000..25ca73c1 --- /dev/null +++ b/libsrc/meshing/surfacegeom.hpp @@ -0,0 +1,70 @@ +#ifndef FILE_SURFACEGEOM +#define FILE_SURFACEGEOM + +/* *************************************************************************/ +/* File: surfacegeom.hpp */ +/* Author: Michael Neunteufel */ +/* Date: Jun. 2020 */ +/* *************************************************************************/ + + +#include + + +namespace netgen +{ + + class DLL_HEADER SurfaceGeometry : public NetgenGeometry + { + function(Point<2>)> func; + double eps=1e-4; + + private: + + void CalcHesse(double u, double v, Vec<3>& f_uu, Vec<3>& f_vv, Vec<3>& f_uv) const; + public: + + SurfaceGeometry(); + SurfaceGeometry(function(Point<2>)> func); + SurfaceGeometry(const SurfaceGeometry& geom); + SurfaceGeometry& operator =(const SurfaceGeometry& geom) + { + func = geom.func; + eps = geom.eps; + return *this; + } + + Array> GetTangentVectors(double u, double v) const; + + virtual Vec<3> GetNormal(int surfind, const Point<3> & p, const PointGeomInfo* gi) const override; + + virtual PointGeomInfo ProjectPoint(int surfind, Point<3> & p) const override; + + virtual void ProjectPointEdge (int surfind, int surfind2, Point<3> & p, + EdgePointGeomInfo* gi = nullptr) const override; + + virtual bool ProjectPointGI (int surfind, Point<3> & p, PointGeomInfo & gi) const override; + + virtual bool CalcPointGeomInfo(int surfind, PointGeomInfo& gi, const Point<3> & p3) const override; + + virtual void PointBetweenEdge(const Point<3> & p1, const Point<3> & p2, double secpoint, + int surfi1, int surfi2, + const EdgePointGeomInfo & ap1, + const EdgePointGeomInfo & ap2, + Point<3> & newp, EdgePointGeomInfo & newgi) const override; + + virtual void PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, + int surfi, + const PointGeomInfo & gi1, + const PointGeomInfo & gi2, + Point<3> & newp, PointGeomInfo & newgi) const override; + + int GenerateMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames); + + }; + +} + + + +#endif //SURFACEGEOM From 8046b19b60aee9c52d2965c90cc82a0699ecd86e Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Thu, 25 Jun 2020 18:39:29 +0200 Subject: [PATCH 17/17] fix facets for 3d bbnd elements --- libsrc/include/nginterface_v2_impl.hpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index 6d6d1b85..20fccd2a 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -111,7 +111,13 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (size_t nr) const ret.faces.num = 0; ret.faces.ptr = NULL; - if (mesh->GetDimension() == 2) + if (mesh->GetDimension() == 3) + { + ret.facets.num = 0; + ret.facets.base = 0; + ret.facets.ptr = nullptr; + } + else if (mesh->GetDimension() == 2) { ret.facets.num = 1; ret.facets.base = 0;