diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 363d6f93..6ecd6266 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1023,6 +1023,11 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) return self.AddFaceDescriptor (fd); }) + .def ("AddSingularity", [](Mesh & self, PointIndex pi, double factor) + { + self[pi].Singularity(factor); + }) + .def ("AddPoints", [](Mesh & self, py::buffer b1) { static Timer timer("Mesh::AddPoints"); @@ -1683,14 +1688,14 @@ project_boundaries : Optional[str] = None }), 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, py::list py_hppts, py::list py_hpbnd) + bool quads, int nx, int ny, bool flip_triangles, py::list py_bbbpts, py::list py_bbbnames, py::list py_hppnts, py::dict/*list*/ py_hpbnd, py::dict py_layers) { 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)); - Array> hppts(py::len(py_hppts)); - Array hpptsfac(py::len(py_hppts)); + Array> hppnts(py::len(py_hppnts)); + Array hppntsfac(py::len(py_hppnts)); Array hpbnd(py::len(py_hpbnd)); Array hpbndfac(py::len(py_hpbnd)); for(int i = 0; i(py::extract(pnt[0])(),py::extract(pnt[1])(),py::extract(pnt[2])()); bbbname[i] = py::extract(py_bbbnames[i])(); } - for(int i = 0; i(py_hppts[i])(); - hppts[i] = Point<3>(py::extract(pnt[0])(),py::extract(pnt[1])(),py::extract(pnt[2])()); - //hpptsfac[i] = py::len(pnt) > 3 ? py::extract(pnt[3])() : 0.0; - hpptsfac[i] = py::extract(pnt[3])(); + py::tuple pnt = py::extract(py_hppnts[i])(); + hppnts[i] = Point<3>(py::extract(pnt[0])(),py::extract(pnt[1])(),py::extract(pnt[2])()); + hppntsfac[i] = py::extract(pnt[3])(); } - for(int i = 0; i(py_hpbnd[i])(); - hpbnd[i] = py::extract(bnd[0])(); - hpbndfac[i] = py::extract(bnd[1])(); - } + int ii=0; + for(auto val : py_hpbnd) + { + hpbnd[ii] = py::cast(val.first); + hpbndfac[ii] = py::cast(val.second); + ii++; + } + + + Array layer_thickness[4]; + bool layer_quad = false; + + for(auto val : py_layers) + { + int index = -1; + if (py::cast(val.first) == "left") index = 0; + else if (py::cast(val.first) == "top") index = 3; + else if (py::cast(val.first) == "right") index = 2; + else if (py::cast(val.first) == "bottom") index = 1; + else if (py::cast(val.first) == "quads") layer_quad = py::cast(val.second); + else throw Exception("Unknown parameter " + string(py::cast(val.first))); + if (index < 0) continue; + + auto list = py::cast(val.second); + layer_thickness[index] = Array(py::len(list)); + for (size_t i = 0; i < py::len(list); i++) + layer_thickness[index][i] = py::cast(list[i]); + } + auto mesh = make_shared(); SetGlobalMesh (mesh); mesh->SetGeometry(geo); ng_geometry = geo; - auto result = geo->GenerateStructuredMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname, hppts, hpptsfac, hpbnd, hpbndfac); + auto result = geo->GenerateStructuredMesh (mesh, quads, nx, ny, flip_triangles, bbbpts, bbbname, hppnts, hppntsfac, hpbnd, hpbndfac, layer_thickness, layer_quad); 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(), py::arg("hppts")=py::list(), py::arg("hpbnd")=py::list()) - ; + }, 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(), py::arg("hppnts")=py::list(), py::arg("hpbnd")=py::dict(), py::arg("boundarylayer")=py::dict());/*, R"raw_string( + Generate a structured 2D surface mesh + + Parameters: + + quads : bool + If True, a quadrilateral mesh is generated. If False, the quads are split to triangles. + + nx : int + Number of cells in x-direction. + + ny : int + Number of cells in y-direction. + + flip_triangles : bool + If set to True together with quads=False the quads are cut the other way round + + bbbpts : list + List of points which should be handled as BBBND and are named with bbbnames. The mesh must be constructed in such a way that the bbbpts coincide with generated points. + + bbbnames : list + List of bbbnd names as strings. Size must coincide with size of bbbpts. + + hppnts : list + If not None it expects a list of the form [ (px1,py1,pz1, hpref1), (px2,py2,pz2, hpref2), ... ] where px,py,pz are the point coordinates which have to be resolved in the mesh and hpref the refinement factor. + + hpbnd : dict + If not None it expects a dictionary of the form {"boundaryname" : hpref } where boundaryname in [left, right, top, bottom] and hpref the refinement factor. + + boundarylayer : dict + If not None it expects a dictionary of the form { "boundaryname" : [t1,...,tn], "quads" : False } where ti denote the thickness of layer i. The number of layers are included in nx/ny. After the layers are placed the remaining number of cells are used to divide the remaining grid uniformly. If quads are set to True quadrilaterals are used inside the boundarylayer. If set False the value of "quads" of the function call is used. + )raw_string");*/ ; py::class_ (m, "ClearSolutionClass") diff --git a/libsrc/meshing/surfacegeom.cpp b/libsrc/meshing/surfacegeom.cpp index d0f3ceb0..a2e81d4d 100644 --- a/libsrc/meshing/surfacegeom.cpp +++ b/libsrc/meshing/surfacegeom.cpp @@ -197,6 +197,7 @@ namespace netgen 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 { newgi.u = ap1.u+secpoint*(ap2.u-ap1.u); @@ -207,6 +208,32 @@ namespace netgen newp = Point<3>(func(Point<2>(newgi.u, newgi.v))); } + + void CheckForBBBPnt(const Array>& bbbpts, const Point<3>& pnt, Array& found, Array& indbbbpts, const Array& pids) + { + 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]; + } + } + } + + void CheckForSingularity(const Array>& hppoints, const Point<3>& pnt, const Array& hpptsfac, shared_ptr & mesh, const Array& pids) + { + for (int k = 0; k < hppoints.Size(); k++) + { + auto diff = pnt - hppoints[k]; + if(diff.Length2() < 1e-14) + { + (*mesh)[pids[pids.Size()-1]].Singularity(hpptsfac[k]); + } + + } + } void SurfaceGeometry :: PointBetween(const Point<3> & p1, const Point<3> & p2, double secpoint, int surfi, @@ -223,7 +250,7 @@ namespace netgen //ProjectPointGI(surfi, newp, newgi); } - int SurfaceGeometry :: GenerateStructuredMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames, const Array>& hppoints, const Array& hpptsfac, const Array& hpbnd, const Array& hpbndfac) + int SurfaceGeometry :: GenerateStructuredMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames, const Array>& hppoints, const Array& hpptsfac, const Array& hpbnd, const Array& hpbndfac, Array layer_thickness[4], bool layer_quad) { mesh->SetDimension(3); @@ -231,41 +258,87 @@ namespace netgen found = false; Array indbbbpts(bbbpts.Size()); + int numx = nx; + int numy = ny; + + size_t total_layer_el[4] = {layer_thickness[0].Size(), layer_thickness[1].Size(), layer_thickness[2].Size(), layer_thickness[3].Size()}; + + double interior_x = 1.0; + double interior_y = 1.0; + for(double scale : layer_thickness[0]) + interior_x -= scale; + for(double scale : layer_thickness[2]) + interior_x -= scale; + for(double scale : layer_thickness[1]) + interior_y -= scale; + for(double scale : layer_thickness[3]) + interior_y -= scale; + + auto AddPoint = [&] (double offsetx, double offsety, Array & pids, Array & pgis) + { + PointGeomInfo pgi; + pgi.trignum = -1; + pgi.u = offsetx; + pgi.v = offsety; + + + Point<3> pnt = Point<3>(func(Point<2>(pgi.u,pgi.v))); + pids.Append(mesh->AddPoint(pnt)); + pgis.Append(pgi); + + CheckForBBBPnt(bbbpts, pnt, found, indbbbpts, pids); + CheckForSingularity(hppoints, pnt, hpptsfac, mesh, pids); + }; + + auto InternalLoop = [&] (double offsety, Array & pids, Array & pgis) + { + int j = 0; + double offsetx = 0.0; + + for(int l=0; l < layer_thickness[0].Size(); l++,j++) + { + AddPoint(offsetx+layer_thickness[0][l]*double(j-l), offsety, pids, pgis); + offsetx += layer_thickness[0][l]; + } + + + for(;j <= nx-total_layer_el[2]; j++) + AddPoint(offsetx + interior_x*double(j-total_layer_el[0])/(nx-total_layer_el[0]-total_layer_el[2]), offsety, pids, pgis); + offsetx += interior_x; + + int startj = j; + for(int l=0; l < layer_thickness[2].Size(); l++, j++) + { + AddPoint(offsetx+layer_thickness[2][layer_thickness[2].Size()-1-l]*double(j-startj-l+1), offsety, pids, pgis); + + offsetx += layer_thickness[2][layer_thickness[2].Size()-1-l]; + } + }; 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; + + int i = 0; + double offsety = 0.0; - 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(int k=0; k < layer_thickness[1].Size(); k++,i++) + { + InternalLoop(offsety, pids, pgis); + offsety += layer_thickness[1][k]; + } - for (int k = 0; k < hppoints.Size(); k++) - { - auto diff = pnt - hppoints[k]; - if(diff.Length2() < 1e-14) - { - (*mesh)[pids[pids.Size()-1]].Singularity(hpptsfac[k]); - } - - } - } + for(; i <= ny-total_layer_el[3]; i++) + { + InternalLoop(offsety, pids, pgis); + offsety += interior_y/(ny-total_layer_el[1]-total_layer_el[3]); + } + offsety -= interior_y/(ny-total_layer_el[1]-total_layer_el[3]); + + for(int k=0; k < layer_thickness[3].Size(); k++,i++) + { + offsety += layer_thickness[3][layer_thickness[3].Size()-1-k]; + InternalLoop(offsety, pids, pgis); + } for (bool f : found) if (!f) @@ -279,14 +352,14 @@ namespace netgen mesh->AddFaceDescriptor(fd); - for(int i=0; i < ny; i++) + for(int i=0; i < numy; i++) { - for(int j=0; j < nx; j++) + for(int j=0; j < numx; j++) { - int base = i * (nx+1) + j; - if (quads) + int base = i * (numx+1) + j; + if (quads || (layer_quad && i < total_layer_el[1]) || (layer_quad && i > numy-1-total_layer_el[3]) || (layer_quad && j < total_layer_el[0]) || (layer_quad && j > numx-1-total_layer_el[2]) ) { - int pnum[4] = {base,base+1,base+nx+2,base+nx+1}; + int pnum[4] = {base,base+1,base+numx+2,base+numx+1}; Element2d el = Element2d(QUAD); for (int i = 0; i < 4; i++) { @@ -305,19 +378,19 @@ namespace netgen { pnum1[0] = base; pnum1[1] = base+1; - pnum1[2] = base+nx+2; + pnum1[2] = base+numx+2; pnum2[0] = base; - pnum2[1] = base+nx+2; - pnum2[2] = base+nx+1; + pnum2[1] = base+numx+2; + pnum2[2] = base+numx+1; } else { pnum1[0] = base; pnum1[1] = base+1; - pnum1[2] = base+nx+1; + pnum1[2] = base+numx+1; pnum2[0] = base+1; - pnum2[1] = base+nx+2; - pnum2[2] = base+nx+1; + pnum2[1] = base+numx+2; + pnum2[2] = base+numx+1; } Element2d el = Element2d(TRIG); @@ -357,7 +430,7 @@ namespace netgen } // needed for codim2 in 3d seg.edgenr = 1; - for(int i=0; i < nx; i++) + for(int i=0; i < numx; i++) { seg[0] = pids[i]; seg[1] = pids[i+1]; @@ -388,18 +461,18 @@ namespace netgen } } - for(int i=0; iAddSegment(seg); @@ -419,18 +492,18 @@ namespace netgen } } - for(int i=0; iAddSegment(seg); @@ -450,18 +523,18 @@ namespace netgen } - for(int i=0; iAddSegment(seg); diff --git a/libsrc/meshing/surfacegeom.hpp b/libsrc/meshing/surfacegeom.hpp index 6402f103..3829be55 100644 --- a/libsrc/meshing/surfacegeom.hpp +++ b/libsrc/meshing/surfacegeom.hpp @@ -62,7 +62,7 @@ namespace netgen const PointGeomInfo & gi2, Point<3> & newp, PointGeomInfo & newgi) const override; - int GenerateStructuredMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames, const Array>& hppoints, const Array& hpptsfac, const Array& hpbnd, const Array& hpbndfac); + int GenerateStructuredMesh(shared_ptr & mesh, bool quads, int nx, int ny, bool flip_triangles, const Array>& bbbpts, const Array& bbbnames, const Array>& hppoints, const Array& hpptsfac, const Array& hpbnd, const Array& hpbndfac, Array layer_thickness[4], bool layer_quad); };