diff --git a/.gitignore b/.gitignore index 36382aa..2b83f83 100644 --- a/.gitignore +++ b/.gitignore @@ -36,3 +36,4 @@ CMakeFiles cmake_install.cmake tmp +.ipynb_checkpoints diff --git a/examples/.gitignore b/examples/.gitignore deleted file mode 100644 index a136337..0000000 --- a/examples/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.pdf diff --git a/examples/jw_mesh_examples.ipynb b/examples/jw_mesh_examples.ipynb index cb7c2f6..7f61537 100644 --- a/examples/jw_mesh_examples.ipynb +++ b/examples/jw_mesh_examples.ipynb @@ -16,18 +16,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 26, "id": "458bfac7-3661-433f-93fc-64d283245089", "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import jw_meshtools as mt\n", - "import meshpy.triangle as triangle\n", - "import numpy.linalg as la\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import numpy.linalg as la\n", "\n", - "length=0.3" + "import meshpy.triangle as triangle\n", + "\n", + "\n", + "length = 0.3" ] }, { @@ -48,21 +50,21 @@ "# Simple mesh rectangle\n", "\n", "# Define closed boundary around a 2D region\n", - "p,v=mt.RectangleSegments([-1.,-1.],[2.5,3.],edge_length=length)\n", - "# Make mesh of this region \n", + "p, v = mt.RectangleSegments([-1.0, -1.0], [2.5, 3.0], edge_length=length)\n", + "# Make mesh of this region\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, edge_length=length)\n", "\n", - "print(\"Points, \",poi[0:5],\"....\",flush=True)\n", - "print(\"Elements, \",tri[0:5],\"...\",flush=True)\n", - "print(\"Boundary edges\",BouE,flush=True)\n", - "print(\"List for boundary edges\",li_BE,flush=True)\n", - "print(\"Boundary triangles\",bou_elem,flush=True)\n", + "print(\"Points, \", poi[0:5], \"....\", flush=True)\n", + "print(\"Elements, \", tri[0:5], \"...\", flush=True)\n", + "print(\"Boundary edges\", BouE, flush=True)\n", + "print(\"List for boundary edges\", li_BE, flush=True)\n", + "print(\"Boundary triangles\", bou_elem, flush=True)\n", "\n", "# Help\n", "print(\"\\n\\n################ Help string :\")\n", - "print( mt.RectangleSegments.__doc__ )\n", - "print( mt.DoTriMesh.__doc__ )" + "print(mt.RectangleSegments.__doc__)\n", + "print(mt.DoTriMesh.__doc__)" ] }, { @@ -73,10 +75,12 @@ "outputs": [], "source": [ "# Simple mesh rectangle with numbers\n", - "p,v=mt.RectangleSegments([-1.,-1.],[2.5,3.],edge_length=3*length)\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=3*length,show=False)\n", + "p, v = mt.RectangleSegments([-1.0, -1.0], [2.5, 3.0], edge_length=3 * length)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(\n", + " p, v, edge_length=3 * length, show=False\n", + ")\n", "\n", - "mt.PlotMeshNumbers(poi,tri)\n" + "mt.PlotMeshNumbers(poi, tri)" ] }, { @@ -94,14 +98,14 @@ "metadata": {}, "outputs": [], "source": [ - "# construct boundary curve from simple lines \n", - "p1,v1=mt.LineSegments([-0.5,0.5],[-1,-1],edge_length=length/5)\n", - "p2,v2=mt.LineSegments([-1,-1],[0.,0.5],edge_length=length/5)\n", - "p3,v3=mt.LineSegments([0.,0.5],[1,1],edge_length=length/7)\n", - "p4,v4=mt.LineSegments([1,1],[-0.5,0.5],edge_length=length/7)\n", - "p,v=mt.AddMultipleSegments(p1,p2,p3,p4)\n", + "# construct boundary curve from simple lines\n", + "p1, v1 = mt.LineSegments([-0.5, 0.5], [-1, -1], edge_length=length / 5)\n", + "p2, v2 = mt.LineSegments([-1, -1], [0.0, 0.5], edge_length=length / 5)\n", + "p3, v3 = mt.LineSegments([0.0, 0.5], [1, 1], edge_length=length / 7)\n", + "p4, v4 = mt.LineSegments([1, 1], [-0.5, 0.5], edge_length=length / 7)\n", + "p, v = mt.AddMultipleSegments(p1, p2, p3, p4)\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length)" + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, edge_length=length)" ] }, { @@ -120,10 +124,9 @@ "outputs": [], "source": [ "# circle as boundary curve\n", - "p,v=mt.CircleSegments([1,2],2,edge_length=length)\n", - "mt.DoTriMesh(p,v,edge_length=length);\n", - "\n", - "print( mt.CircleSegments.__doc__ )" + "p, v = mt.CircleSegments([1, 2], 2, edge_length=length)\n", + "mt.DoTriMesh(p, v, edge_length=length)\n", + "print(mt.CircleSegments.__doc__)" ] }, { @@ -142,10 +145,10 @@ "outputs": [], "source": [ "#\n", - "p1,v1=mt.LineSegments([2,2],[-1,-3],edge_length=length)\n", - "p2,v2=mt.LineSegments([-1,-4],[3,-1],num_points=10)\n", - "p,v=mt.AddSegments(p1,p2,closed=True)\n", - "mt.DoTriMesh(p,v,edge_length=length);" + "p1, v1 = mt.LineSegments([2, 2], [-1, -3], edge_length=length)\n", + "p2, v2 = mt.LineSegments([-1, -4], [3, -1], num_points=10)\n", + "p, v = mt.AddSegments(p1, p2, closed=True)\n", + "mt.DoTriMesh(p, v, edge_length=length);" ] }, { @@ -166,8 +169,8 @@ "#\n", "# rectangle with smooth corners\n", "#\n", - "p,v=mt.ORecSegments([1,2],[7,6],0.3,edge_length=length,num_pc=10)\n", - "mt.DoTriMesh(p,v,edge_length=length);" + "p, v = mt.ORecSegments([1, 2], [7, 6], 0.3, edge_length=length, num_pc=10)\n", + "mt.DoTriMesh(p, v, edge_length=length);" ] }, { @@ -185,14 +188,18 @@ "metadata": {}, "outputs": [], "source": [ - "# \n", + "#\n", "# two semicircles\n", "#\n", - "p1,v1=mt.CircleSegments([1.,0],1,a_min=-np.pi/2,a_max=np.pi/2,num_points=20)\n", - "p2,v2=mt.CircleSegments([1,0],3,a_min=np.pi/2.,a_max=3.*np.pi/2,num_points=20)\n", - "p,v=mt.AddSegments(p1,p2,closed=True)\n", - "# plot mesh \n", - "mt.DoTriMesh(p,v,edge_length=length);" + "p1, v1 = mt.CircleSegments(\n", + " [1.0, 0], 1, a_min=-np.pi / 2, a_max=np.pi / 2, num_points=20\n", + ")\n", + "p2, v2 = mt.CircleSegments(\n", + " [1, 0], 3, a_min=np.pi / 2.0, a_max=3.0 * np.pi / 2, num_points=20\n", + ")\n", + "p, v = mt.AddSegments(p1, p2, closed=True)\n", + "# plot mesh\n", + "mt.DoTriMesh(p, v, edge_length=length);" ] }, { @@ -213,13 +220,13 @@ "#\n", "# boundary curve defined by simple points\n", "#\n", - "t=np.linspace(0,2*np.pi,120)\n", - "r=3+np.sin(8*t)\n", - "x=r*np.cos(t)\n", - "y=r*np.sin(t)\n", - "p=[(x[j],y[j]) for j in range(len(t))]\n", - "p1,v1=mt.PointSegments(p)\n", - "mt.DoTriMesh(p1,v1,edge_length=length);" + "t = np.linspace(0, 2 * np.pi, 120)\n", + "r = 3 + np.sin(8 * t)\n", + "x = r * np.cos(t)\n", + "y = r * np.sin(t)\n", + "p = [(x[j], y[j]) for j in range(len(t))]\n", + "p1, v1 = mt.PointSegments(p)\n", + "mt.DoTriMesh(p1, v1, edge_length=length);" ] }, { @@ -241,37 +248,42 @@ "# Example for using directly triangle\n", "#\n", "\n", + "\n", "def round_trip_connect(start, end):\n", - " return [(i, i+1) for i in range(start, end)] + [(end, start)]\n", + " return [(i, i + 1) for i in range(start, end)] + [(end, start)]\n", "\n", - "points = [(1,0),(1,1),(-1,1),(-1,-1),(1,-1),(1,0)]\n", - "facets = round_trip_connect(0, len(points)-1)\n", + "\n", + "points = [(1, 0), (1, 1), (-1, 1), (-1, -1), (1, -1), (1, 0)]\n", + "facets = round_trip_connect(0, len(points) - 1)\n", "\n", "circ_start = len(points)\n", "points.extend(\n", - " (3 * np.cos(angle), 3 * np.sin(angle))\n", - " for angle in np.linspace(0, 2*np.pi, 29, endpoint=False))\n", + " (3 * np.cos(angle), 3 * np.sin(angle))\n", + " for angle in np.linspace(0, 2 * np.pi, 29, endpoint=False)\n", + ")\n", + "\n", + "facets.extend(round_trip_connect(circ_start, len(points) - 1))\n", "\n", - "facets.extend(round_trip_connect(circ_start, len(points)-1))\n", "\n", "def needs_refinement(vertices, area):\n", - " bary = np.sum(np.array(vertices), axis=0)/3\n", - " max_area = 0.01 + abs((la.norm(bary, np.inf)-1))*0.1\n", + " bary = np.sum(np.array(vertices), axis=0) / 3\n", + " max_area = 0.01 + abs(la.norm(bary, np.inf) - 1) * 0.1\n", " return bool(area > max_area)\n", "\n", + "\n", "info = triangle.MeshInfo()\n", "info.set_points(points)\n", - "info.set_holes([(0,0)])\n", + "info.set_holes([(0, 0)])\n", "info.set_facets(facets)\n", "\n", "mesh = triangle.build(info, refinement_func=needs_refinement)\n", - "#mesh = triangle.build(info) \n", + "# mesh = triangle.build(info)\n", "\n", "mesh_points = np.array(mesh.points)\n", "mesh_tris = np.array(mesh.elements)\n", "\n", - "print(mesh_points[0:5],\"...\")\n", - "print(mesh_tris[0:5],\"....\")\n", + "print(mesh_points[0:5], \"...\")\n", + "print(mesh_tris[0:5], \"....\")\n", "plt.triplot(mesh_points[:, 0], mesh_points[:, 1], mesh_tris)\n", "plt.show()" ] @@ -294,18 +306,17 @@ "#\n", "# rectangle and inner circle\n", "#\n", - "p1,v1=mt.RectangleSegments([-2,-2],[2.5,3],edge_length=length)\n", + "p1, v1 = mt.RectangleSegments([-2, -2], [2.5, 3], edge_length=length)\n", "\n", - "p2,v2=mt.CircleSegments([1,1],1,edge_length=length/5)\n", - "p,v=mt.AddCurves(p1,v1,p2,v2)\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length)\n", - "print(\"Points, \",poi[0:5],\"...\")\n", - "print(\"Elements, \",tri[0:5],\"...\")\n", - "print(\"Boundary Edges\",BouE[0:5],\"...\")\n", - "print(\"List boundary edges\",li_BE)\n", - "print(\"Inner Curves\",CuE[0:5],\"...\")\n", - "print(\"List inner Curve\",li_CE)\n", - "\n" + "p2, v2 = mt.CircleSegments([1, 1], 1, edge_length=length / 5)\n", + "p, v = mt.AddCurves(p1, v1, p2, v2)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, edge_length=length)\n", + "print(\"Points, \", poi[0:5], \"...\")\n", + "print(\"Elements, \", tri[0:5], \"...\")\n", + "print(\"Boundary Edges\", BouE[0:5], \"...\")\n", + "print(\"List boundary edges\", li_BE)\n", + "print(\"Inner Curves\", CuE[0:5], \"...\")\n", + "print(\"List inner Curve\", li_CE)" ] }, { @@ -326,18 +337,18 @@ "#\n", "# rectangle and inner line\n", "#\n", - "p1,v1=mt.RectangleSegments([-2,-2],[2.5,3],edge_length=length)\n", - "p2,v2=mt.LineSegments([0,0],[1,1],edge_length=length/5)\n", + "p1, v1 = mt.RectangleSegments([-2, -2], [2.5, 3], edge_length=length)\n", + "p2, v2 = mt.LineSegments([0, 0], [1, 1], edge_length=length / 5)\n", "\n", "\n", - "p3,v3=mt.LineSegments([-1,1],[0,-1],edge_length=length/5)\n", - "p4,v4=mt.LineSegments([0,-1],[1,-1],edge_length=length/5)\n", + "p3, v3 = mt.LineSegments([-1, 1], [0, -1], edge_length=length / 5)\n", + "p4, v4 = mt.LineSegments([0, -1], [1, -1], edge_length=length / 5)\n", "# connect line 3 and 4 first\n", - "p5,v5=mt.AddSegments(p3,p4)\n", + "p5, v5 = mt.AddSegments(p3, p4)\n", "\n", - "p,v,indizes=mt.AddMultipleCurves(p1,v1,p2,v2,p5,v5)\n", + "p, v, indizes = mt.AddMultipleCurves(p1, v1, p2, v2, p5, v5)\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length)" + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, edge_length=length)" ] }, { @@ -356,27 +367,27 @@ "outputs": [], "source": [ "#\n", - "# rectangle and more complicated inner curves \n", + "# rectangle and more complicated inner curves\n", "#\n", - "p1,v1=mt.RectangleSegments([-2,-2],[2.5,3],edge_length=length)\n", + "p1, v1 = mt.RectangleSegments([-2, -2], [2.5, 3], edge_length=length)\n", "\n", - "p2,v2=mt.CircleSegments([1,1],1,edge_length=length/5)\n", - "p,v=mt.AddCurves(p1,v1,p2,v2)\n", + "p2, v2 = mt.CircleSegments([1, 1], 1, edge_length=length / 5)\n", + "p, v = mt.AddCurves(p1, v1, p2, v2)\n", "\n", "# use connect if segments might have nearly the same points\n", - "p3,v3=mt.LineSegments([-1,-2],[-1,3],edge_length=length/4)\n", - "p,v=mt.AddCurves(p,v,p3,v3,connect=True,eps=1e-12)\n", + "p3, v3 = mt.LineSegments([-1, -2], [-1, 3], edge_length=length / 4)\n", + "p, v = mt.AddCurves(p, v, p3, v3, connect=True, eps=1e-12)\n", "\n", - "p4,v4=mt.LineSegments([-1,1],[0,1],edge_length=length/5)\n", - "p,v=mt.AddCurves(p,v,p4,v4,connect=True,eps=1e-12)\n", + "p4, v4 = mt.LineSegments([-1, 1], [0, 1], edge_length=length / 5)\n", + "p, v = mt.AddCurves(p, v, p4, v4, connect=True, eps=1e-12)\n", "\n", "# or shift inner curve slightly\n", "\n", - "epsilon=1e-6\n", - "p5,v5=mt.LineSegments([1,-2+epsilon],[3-epsilon,-1],edge_length=length/5)\n", - "p,v=mt.AddCurves(p,v,p5,v5)\n", + "epsilon = 1e-6\n", + "p5, v5 = mt.LineSegments([1, -2 + epsilon], [3 - epsilon, -1], edge_length=length / 5)\n", + "p, v = mt.AddCurves(p, v, p5, v5)\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length)\n" + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, edge_length=length)" ] }, { @@ -396,18 +407,20 @@ "source": [ "#\n", "# rectangle with holes\n", - "p1,v1=mt.LineSegments([-2,-3],[2,-3],num_points=12)\n", - "p2,v2=mt.LineSegments([2,3],[-2,3],num_points=12)\n", - "p,v=mt.AddSegments(p1,p2,closed=True)\n", + "p1, v1 = mt.LineSegments([-2, -3], [2, -3], num_points=12)\n", + "p2, v2 = mt.LineSegments([2, 3], [-2, 3], num_points=12)\n", + "p, v = mt.AddSegments(p1, p2, closed=True)\n", "\n", "# define the boundary curves of holes\n", - "p3,v3=mt.CircleSegments([-0.5,0.5],0.5,edge_length=length)\n", - "p,v=mt.AddCurves(p,v,p3,v3)\n", - "p4,v4=mt.CircleSegments([1,-1],0.5,edge_length=length)\n", - "p,v=mt.AddCurves(p,v,p4,v4)\n", + "p3, v3 = mt.CircleSegments([-0.5, 0.5], 0.5, edge_length=length)\n", + "p, v = mt.AddCurves(p, v, p3, v3)\n", + "p4, v4 = mt.CircleSegments([1, -1], 0.5, edge_length=length)\n", + "p, v = mt.AddCurves(p, v, p4, v4)\n", "\n", - "# the array holes contain points in the regions to be removed \n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length,holes=[(-0.4,0.4),(0.95,-0.8)])" + "# the array holes contain points in the regions to be removed\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(\n", + " p, v, edge_length=length, holes=[(-0.4, 0.4), (0.95, -0.8)]\n", + ")" ] }, { @@ -427,29 +440,31 @@ "source": [ "# boundary nodes\n", "# rectangle with holes\n", - "p,v=mt.RectangleSegments([-2,-3],[2,3],edge_length=length)\n", + "p, v = mt.RectangleSegments([-2, -3], [2, 3], edge_length=length)\n", "\n", - "p3,v3=mt.CircleSegments([-0.5,0.5],0.5,edge_length=length)\n", - "p4,v4=mt.CircleSegments([1,-1],0.5,edge_length=length)\n", - "p,v,ii=mt.AddMultipleCurves(p,v,p3,v3,p4,v4)\n", + "p3, v3 = mt.CircleSegments([-0.5, 0.5], 0.5, edge_length=length)\n", + "p4, v4 = mt.CircleSegments([1, -1], 0.5, edge_length=length)\n", + "p, v, ii = mt.AddMultipleCurves(p, v, p3, v3, p4, v4)\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length,holes=[(-0.4,0.4),(0.95,-0.8)],show=False)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(\n", + " p, v, edge_length=length, holes=[(-0.4, 0.4), (0.95, -0.8)], show=False\n", + ")\n", "\n", "# node numbers used for search\n", "all_nodes = np.arange(len(poi))\n", "# points to be found\n", - "p0=[ [2*np.cos(t),2*np.sin(t)] for t in np.linspace(0,2*np.pi,15)]\n", + "p0 = [[2 * np.cos(t), 2 * np.sin(t)] for t in np.linspace(0, 2 * np.pi, 15)]\n", "\n", "# search\n", - "nn,dd = mt.FindClosestNode(all_nodes,poi,p0)\n", + "nn, dd = mt.FindClosestNode(all_nodes, poi, p0)\n", "\n", "\n", - "print(\"Points given\\n\",p0)\n", - "print(\"Node number\\n\",nn)\n", - "print(\"Distance from p0\\n\",dd)\n", + "print(\"Points given\\n\", p0)\n", + "print(\"Node number\\n\", nn)\n", + "print(\"Distance from p0\\n\", dd)\n", "\n", "plt.triplot(poi[:, 0], poi[:, 1], tri)\n", - "plt.plot(poi[nn,0],poi[nn,1],'or');" + "plt.plot(poi[nn, 0], poi[nn, 1], \"or\");" ] }, { @@ -469,15 +484,15 @@ "source": [ "# take mesh from above\n", "\n", - "# Find boundary nodes/segments between the points below \n", + "# Find boundary nodes/segments between the points below\n", "# Two types of boundaries, Nodes or segments\n", - "# Ps: \n", - "Ps = [ [-2,3],[2,-3],[-2,3],[-0.5,0.5],[-0.5,0.5],[1,-1],[1,-1] ]\n", - "Ps_types = ['Nodes','Segments','Nodes','Segments']\n", + "# Ps:\n", + "Ps = [[-2, 3], [2, -3], [-2, 3], [-0.5, 0.5], [-0.5, 0.5], [1, -1], [1, -1]]\n", + "Ps_types = [\"Nodes\", \"Segments\", \"Nodes\", \"Segments\"]\n", "\n", - "bseg=mt.RetrieveSegments(poi,BouE,li_BE,Ps,Ps_types)\n", - "for i in range( len(Ps_types) ):\n", - " print(\"bseg[\",i,\"] : \",Ps_types[i],\"\\n\",bseg[i])\n", + "bseg = mt.RetrieveSegments(poi, BouE, li_BE, Ps, Ps_types)\n", + "for i in range(len(Ps_types)):\n", + " print(\"bseg[\", i, \"] : \", Ps_types[i], \"\\n\", bseg[i])\n", "\n", "# !!!!!!!!!!!!!!!!!!!!!!!!!\n", "# bseg[0] contains all nodes (Ps_types[0]) between Ps[0] and Ps[1]\n", @@ -490,7 +505,7 @@ "plt.triplot(poi[:, 0], poi[:, 1], tri)\n", "\n", "for i in range(len(Ps_types)):\n", - " mt.PlotBoundary(poi,bseg[i],Ps_types[i])\n", + " mt.PlotBoundary(poi, bseg[i], Ps_types[i])\n", "plt.show()" ] }, @@ -511,24 +526,26 @@ "source": [ "#\n", "# rectangle with holes\n", - "p,v=mt.RectangleSegments([-2.,-3.],[2,3.],edge_length=length)\n", - "p3,v3=mt.CircleSegments([-0.5,0.5],0.5,edge_length=length/4)\n", - "p,v=mt.AddCurves(p,v,p3,v3)\n", - "p4,v4=mt.CircleSegments([1,-1],0.5,edge_length=length/4)\n", - "p,v=mt.AddCurves(p,v,p4,v4)\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length,show=False)\n", + "p, v = mt.RectangleSegments([-2.0, -3.0], [2, 3.0], edge_length=length)\n", + "p3, v3 = mt.CircleSegments([-0.5, 0.5], 0.5, edge_length=length / 4)\n", + "p, v = mt.AddCurves(p, v, p3, v3)\n", + "p4, v4 = mt.CircleSegments([1, -1], 0.5, edge_length=length / 4)\n", + "p, v = mt.AddCurves(p, v, p4, v4)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(\n", + " p, v, edge_length=length, show=False\n", + ")\n", "\n", "# same as before, with CuE and li_DE\n", - "Ps = [[-0.5,0.5],[-0.5,0.5],[1,-1],[1,-1] ]\n", - "Ps_types = ['Nodes','Segments']\n", + "Ps = [[-0.5, 0.5], [-0.5, 0.5], [1, -1], [1, -1]]\n", + "Ps_types = [\"Nodes\", \"Segments\"]\n", "\n", - "cseg=mt.RetrieveSegments(poi,CuE,li_CE,Ps,Ps_types)\n", - "print(\"cseg\",cseg)\n", + "cseg = mt.RetrieveSegments(poi, CuE, li_CE, Ps, Ps_types)\n", + "print(\"cseg\", cseg)\n", "\n", "plt.triplot(poi[:, 0], poi[:, 1], tri)\n", "\n", - "mt.PlotBoundary(poi,cseg[0],'Nodes')\n", - "mt.PlotBoundary(poi,cseg[1],'Segments')\n", + "mt.PlotBoundary(poi, cseg[0], \"Nodes\")\n", + "mt.PlotBoundary(poi, cseg[1], \"Segments\")\n", "plt.show()" ] }, @@ -550,28 +567,28 @@ "#\n", "# rectangle and inner line\n", "#\n", - "p1,v1=mt.RectangleSegments([-2,-2],[2.5,3],edge_length=length)\n", + "p1, v1 = mt.RectangleSegments([-2, -2], [2.5, 3], edge_length=length)\n", "\n", - "p2,v2=mt.CircleSegments([1,1],1,edge_length=length/5)\n", - "p,v=mt.AddCurves(p1,v1,p2,v2)\n", + "p2, v2 = mt.CircleSegments([1, 1], 1, edge_length=length / 5)\n", + "p, v = mt.AddCurves(p1, v1, p2, v2)\n", "\n", "# use connect if segments might have nearly the same points\n", - "p3,v3=mt.LineSegments([-1,-2],[-1,3],edge_length=length/4)\n", - "p,v=mt.AddCurves(p,v,p3,v3,connect=True,eps=1e-12)\n", + "p3, v3 = mt.LineSegments([-1, -2], [-1, 3], edge_length=length / 4)\n", + "p, v = mt.AddCurves(p, v, p3, v3, connect=True, eps=1e-12)\n", "\n", - "p4,v4=mt.LineSegments([-1,1],[0,1],edge_length=length/5)\n", - "p,v=mt.AddCurves(p,v,p4,v4,connect=True,eps=1e-12)\n", + "p4, v4 = mt.LineSegments([-1, 1], [0, 1], edge_length=length / 5)\n", + "p, v = mt.AddCurves(p, v, p4, v4, connect=True, eps=1e-12)\n", "\n", "# or shift inner curve slightly away from existing points/curves\n", - "epsilon=1e-6\n", - "p5,v5=mt.LineSegments([1,-2+epsilon],[3-epsilon,-1],edge_length=length/5)\n", - "p,v=mt.AddCurves(p,v,p5,v5)\n", + "epsilon = 1e-6\n", + "p5, v5 = mt.LineSegments([1, -2 + epsilon], [3 - epsilon, -1], edge_length=length / 5)\n", + "p, v = mt.AddCurves(p, v, p5, v5)\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, edge_length=length)\n", "\n", "# plot all boundaries and inner curves\n", - "mt.PlotBoundary(poi,BouE,'Segments')\n", - "mt.PlotBoundary(poi,CuE,'Curves')" + "mt.PlotBoundary(poi, BouE, \"Segments\")\n", + "mt.PlotBoundary(poi, CuE, \"Curves\")" ] }, { @@ -590,14 +607,14 @@ "outputs": [], "source": [ "#\n", - "# rectangle and local refinement \n", + "# rectangle and local refinement\n", "#\n", - "p1,v1=mt.RectangleSegments([0,0],[1,1],num_points=100)\n", - "p2,v2=mt.RectangleSegments([0.05,0.05],[0.95,0.95],num_points=40)\n", - "p,v=mt.AddCurves(p1,v1,p2,v2)\n", - "p3,v3=mt.RectangleSegments([0.1,0.1],[0.9,0.9],num_points=20)\n", - "p,v=mt.AddCurves(p,v,p3,v3)\n", - "mt.DoTriMesh(p,v,edge_length=length);" + "p1, v1 = mt.RectangleSegments([0, 0], [1, 1], num_points=100)\n", + "p2, v2 = mt.RectangleSegments([0.05, 0.05], [0.95, 0.95], num_points=40)\n", + "p, v = mt.AddCurves(p1, v1, p2, v2)\n", + "p3, v3 = mt.RectangleSegments([0.1, 0.1], [0.9, 0.9], num_points=20)\n", + "p, v = mt.AddCurves(p, v, p3, v3)\n", + "mt.DoTriMesh(p, v, edge_length=length);" ] }, { @@ -618,26 +635,28 @@ "#\n", "# 2D curve with local mesh refinement I\n", "#\n", - "# \n", - "t=np.linspace(0,2*np.pi,120)\n", - "r=3+np.sin(8*t)\n", - "x=r*np.cos(t)\n", - "y=r*np.sin(t)\n", - "p=[(x[j],y[j]) for j in range(len(t))]\n", - "p1,v1=mt.PointSegments(p)\n", + "#\n", + "t = np.linspace(0, 2 * np.pi, 120)\n", + "r = 3 + np.sin(8 * t)\n", + "x = r * np.cos(t)\n", + "y = r * np.sin(t)\n", + "p = [(x[j], y[j]) for j in range(len(t))]\n", + "p1, v1 = mt.PointSegments(p)\n", "# function for refinement\n", "\n", + "\n", "def myrefine1(tri_points, area):\n", - " center_tri = np.sum(np.array(tri_points), axis=0)/3.\n", - " x=center_tri[0]\n", - " y=center_tri[1]\n", - " if x>0:\n", - " max_area=0.05*(1+3*x)\n", - " else:\n", - " max_area=0.05\n", - " return bool(area>max_area)\n", + " center_tri = np.sum(np.array(tri_points), axis=0) / 3.0\n", + " x = center_tri[0]\n", + " _y = center_tri[1]\n", + " if x > 0:\n", + " max_area = 0.05 * (1 + 3 * x)\n", + " else:\n", + " max_area = 0.05\n", + " return bool(area > max_area)\n", + "\n", "\n", - "mt.DoTriMesh(p1,v1,tri_refine=myrefine1);" + "mt.DoTriMesh(p1, v1, tri_refine=myrefine1);" ] }, { @@ -649,11 +668,13 @@ "source": [ "# function for refinement\n", "def myrefine2(tri_points, area):\n", - " center_tri = np.sum(np.array(tri_points), axis=0)/3.\n", - " r=np.sqrt(center_tri[0]**2+center_tri[1]**2) \n", - " max_area=0.3+(0.01-0.3)/(1+0.5*(r-1)**2)\n", - " return bool(area>max_area);\n", - "mt.DoTriMesh(p1,v1,tri_refine=myrefine2);" + " center_tri = np.sum(np.array(tri_points), axis=0) / 3.0\n", + " r = np.sqrt(center_tri[0] ** 2 + center_tri[1] ** 2)\n", + " max_area = 0.3 + (0.01 - 0.3) / (1 + 0.5 * (r - 1) ** 2)\n", + " return bool(area > max_area)\n", + "\n", + "\n", + "mt.DoTriMesh(p1, v1, tri_refine=myrefine2);" ] }, { @@ -676,17 +697,18 @@ "# !! 2 plots\n", "#\n", "# take p1 from above\n", - "p2,v2=mt.CircleSegments([0,0],1,edge_length=0.05)\n", - "p,v=mt.AddCurves(p1,v1,p2,v2)\n", + "p2, v2 = mt.CircleSegments([0, 0], 1, edge_length=0.05)\n", + "p, v = mt.AddCurves(p1, v1, p2, v2)\n", + "\n", "\n", "def myrefine3(tri_points, area):\n", - " center_tri = np.sum(np.array(tri_points), axis=0)/3.\n", - " r=np.sqrt(center_tri[0]**2+center_tri[1]**2) \n", - " max_area=0.4+(0.01-0.3)/(1+0.5*(r-1)**2)\n", - " return bool(area>max_area);\n", - " \n", + " center_tri = np.sum(np.array(tri_points), axis=0) / 3.0\n", + " r = np.sqrt(center_tri[0] ** 2 + center_tri[1] ** 2)\n", + " max_area = 0.4 + (0.01 - 0.3) / (1 + 0.5 * (r - 1) ** 2)\n", + " return bool(area > max_area)\n", + "\n", "\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,tri_refine=myrefine3)" + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(p, v, tri_refine=myrefine3)" ] }, { @@ -704,26 +726,30 @@ "metadata": {}, "outputs": [], "source": [ + "from scipy.spatial import cKDTree\n", + "\n", + "\n", "#\n", "# 2D curve with local refinement III\n", - "# \n", + "#\n", "#\n", "# take p1 from above\n", - "nodes=range(len(p1))\n", + "nodes = range(len(p1))\n", "# define tree to speed up node search\n", - "from scipy.spatial import cKDTree\n", - "p1tree=cKDTree(np.array(p1))\n", + "p1tree = cKDTree(np.array(p1))\n", + "\n", "\n", "# function for refinement\n", "def myrefine3(tri_points, area):\n", - " center_tri = np.sum(np.array(tri_points), axis=0)/3.\n", - " p0=[(center_tri[0],center_tri[1])]\n", - " node,r=mt.FindClosestNode(nodes,[],p0,tree=p1tree)\n", - " r=r[0]\n", - " max_area=0.3+(0.01-0.3)/(1+r**2) \n", - " return bool(area>max_area);\n", + " center_tri = np.sum(np.array(tri_points), axis=0) / 3.0\n", + " p0 = [(center_tri[0], center_tri[1])]\n", + " _node, r = mt.FindClosestNode(nodes, [], p0, tree=p1tree)\n", + " r = r[0]\n", + " max_area = 0.3 + (0.01 - 0.3) / (1 + r**2)\n", + " return bool(area > max_area)\n", "\n", - "mt.DoTriMesh(p1,v1,tri_refine=myrefine3);" + "\n", + "mt.DoTriMesh(p1, v1, tri_refine=myrefine3);" ] }, { @@ -742,17 +768,19 @@ "outputs": [], "source": [ "# Simple mesh rectangle with second order points\n", - "p,v=mt.RectangleSegments([-1.,-1.],[2.5,3.],edge_length=length)\n", - "poi,tri,BouE,li_BE,bou_elem,CuE,li_CE=mt.DoTriMesh(p,v,edge_length=length,order=2,show=None)\n", - "\n", - "plt.triplot(poi[:, 0], poi[:, 1], tri[:,0:3])\n", - "maxi=np.max(tri[:,0:3])+1\n", - "plt.plot(poi[maxi:,0],poi[maxi:,1],'k*')\n", - "mt.PlotBoundary(poi,np.array(BouE),'Segments') \n", + "p, v = mt.RectangleSegments([-1.0, -1.0], [2.5, 3.0], edge_length=length)\n", + "poi, tri, BouE, li_BE, bou_elem, CuE, li_CE = mt.DoTriMesh(\n", + " p, v, edge_length=length, order=2, show=None\n", + ")\n", + "\n", + "plt.triplot(poi[:, 0], poi[:, 1], tri[:, 0:3])\n", + "maxi = np.max(tri[:, 0:3]) + 1\n", + "plt.plot(poi[maxi:, 0], poi[maxi:, 1], \"k*\")\n", + "mt.PlotBoundary(poi, np.array(BouE), \"Segments\")\n", "plt.show()\n", - "print(\"points:\",poi[0:5],\"....\")\n", - "print(\"elements\",tri[0:5],\"....\")\n", - "print(\"boundary\",BouE)" + "print(\"points:\", poi[0:5], \"....\")\n", + "print(\"elements\", tri[0:5], \"....\")\n", + "print(\"boundary\", BouE)" ] }, { @@ -773,54 +801,40 @@ "# connect mesh\n", "\n", "# mesh A\n", - "p1,v1=mt.LineSegments([0,1],[0,0],edge_length=length)\n", - "p2,v2=mt.LineSegments([0,0],[1,0],edge_length=length)\n", - "p,v=mt.AddSegments(p1,p2)\n", - "p1,v1=mt.CircleSegments([0,0],1,a_min=0,a_max=np.pi/2,edge_length=length)\n", - "p,v=mt.AddSegments(p,p1)\n", - "pA,tA,bA,lA,bou_elemA,cuA,lcA=mt.DoTriMesh(p,v,edge_length=length)\n", - "\n", - "#mesh B\n", - "p1,v1=mt.CircleSegments([0,0],1,a_min=0,a_max=np.pi/2,edge_length=length)\n", - "p2,v2=mt.LineSegments([0,1],[2,2],edge_length=length)\n", - "p,v=mt.AddSegments(p1,p2)\n", - "p1,v1=mt.CircleSegments([0,0],2*np.sqrt(2),a_min=np.pi/4,a_max=0,edge_length=length)\n", - "p,v=mt.AddSegments(p,p1)\n", - "p1,v1=mt.LineSegments([2*np.sqrt(2),0],[1,0],edge_length=length)\n", - "p,v=mt.AddSegments(p,p1)\n", - "pB,tB,bB,lB,bou_elemB,cuB,lcB=mt.DoTriMesh(p,v,edge_length=length)\n", - "\n", - "#connect\n", - "p,t,b,bl,idn=mt.ConnectMesh(pA,tA,bA,pB,tB,bB,epsilon=1e-8)\n", - "plt.triplot(p[:, 0],p[:, 1],t[:,0:3])\n", - "k=[x[0] for x in idn]\n", - "plt.plot(p[k,0],p[k,1],'ro',mfc='none')\n", - "\n", - "mt.PlotBoundary(p,np.array(b),'Segments')\n", + "p1, v1 = mt.LineSegments([0, 1], [0, 0], edge_length=length)\n", + "p2, v2 = mt.LineSegments([0, 0], [1, 0], edge_length=length)\n", + "p, v = mt.AddSegments(p1, p2)\n", + "p1, v1 = mt.CircleSegments([0, 0], 1, a_min=0, a_max=np.pi / 2, edge_length=length)\n", + "p, v = mt.AddSegments(p, p1)\n", + "pA, tA, bA, lA, bou_elemA, cuA, lcA = mt.DoTriMesh(p, v, edge_length=length)\n", + "\n", + "# mesh B\n", + "p1, v1 = mt.CircleSegments([0, 0], 1, a_min=0, a_max=np.pi / 2, edge_length=length)\n", + "p2, v2 = mt.LineSegments([0, 1], [2, 2], edge_length=length)\n", + "p, v = mt.AddSegments(p1, p2)\n", + "p1, v1 = mt.CircleSegments(\n", + " [0, 0], 2 * np.sqrt(2), a_min=np.pi / 4, a_max=0, edge_length=length\n", + ")\n", + "p, v = mt.AddSegments(p, p1)\n", + "p1, v1 = mt.LineSegments([2 * np.sqrt(2), 0], [1, 0], edge_length=length)\n", + "p, v = mt.AddSegments(p, p1)\n", + "pB, tB, bB, lB, bou_elemB, cuB, lcB = mt.DoTriMesh(p, v, edge_length=length)\n", + "\n", + "# connect\n", + "p, t, b, bl, idn = mt.ConnectMesh(pA, tA, bA, pB, tB, bB, epsilon=1e-8)\n", + "plt.triplot(p[:, 0], p[:, 1], t[:, 0:3])\n", + "k = [x[0] for x in idn]\n", + "plt.plot(p[k, 0], p[k, 1], \"ro\", mfc=\"none\")\n", + "\n", + "mt.PlotBoundary(p, np.array(b), \"Segments\")\n", "plt.show()\n", "\n", - "Ps=[ [1,0],[1,0] ]\n", - "bseg=mt.RetrieveSegments(p,b,bl,Ps,['Nodes'])\n", - "mt.PlotBoundary(p,bseg[0],'Nodes')\n", + "Ps = [[1, 0], [1, 0]]\n", + "bseg = mt.RetrieveSegments(p, b, bl, Ps, [\"Nodes\"])\n", + "mt.PlotBoundary(p, bseg[0], \"Nodes\")\n", "\n", - "plt.show() " + "plt.show()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "02e4c20b-52e1-44a3-93fe-f598a8be6cbb", - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c6cf20a4-4895-4ef3-9509-03997eb60262", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -839,7 +853,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.2" + "version": "3.12.4" } }, "nbformat": 4, diff --git a/examples/jw_meshtools.py b/examples/jw_meshtools.py index e023f8b..c9ed54d 100644 --- a/examples/jw_meshtools.py +++ b/examples/jw_meshtools.py @@ -1754,7 +1754,7 @@ def AddCurves(p1, v1, p2, v2, connect=False, connect_points=None, eps=1e-12): """ p,v = AddCurves(p1,v1,p2,v2,connect=False,connect_points=[],eps=1e-12) """ - if connect_points: + if connect_points is None: connect_points = [] # make one list diff --git a/pyproject.toml b/pyproject.toml index d162bf0..f6f5671 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -79,14 +79,11 @@ extend-ignore = [ "E242", # tab after comma "E265", # block comment should start with # "E402", # module level import not at the top of file - - # TODO - "UP031", # f-string ] [tool.ruff.lint.per-file-ignores] -"examples/jw_meshtools.py" = ["N806", "N803", "N802"] -"examples/jw_mesh_examples.ipynb" = ["E", "Q", "N", "I", "UP", "W", "F"] +"examples/jw_meshtools.py" = ["N802", "N803", "N806"] +"examples/jw_mesh_examples.ipynb" = ["N802", "N803", "N806", "N816"] [tool.ruff.lint.flake8-quotes] inline-quotes = "double"