update

app/HVSCMesh/area_adaptive/rreo_EPI.py

@@ -2,8 +2,10 @@
 import gmsh
 import matplotlib.pyplot as plt
 
+#from TriRadiusRatio import TriRadiusRatio
+from fealpy.experimental.mesh.mesh_quality import RadiusRatioQuality
+
 from fealpy.mesh import TriangleMesh
-from TriRadiusRatio import TriRadiusRatio
 
 mesh = TriangleMesh.from_meshio('area_EPI_adaptive.vtu')
 node = mesh.entity('node')
@@ -15,21 +17,31 @@
 max_angle = np.max(angle,axis=1)
 angles = max_angle*(180/np.pi)
 m90 = np.sum(angles>90)
+mesh.celldata['angles'] = angles
 
 fig,axes1= plt.subplots()
 axes1.set_ylim(0,8500)
 fig.text(0.15,0.85,f'Number of angles>90:{m90}',fontsize=12)
 mesh.show_angle(axes1,max_angle)
 plt.show()
 
-mesh.celldata['angles'] = angles
-mesh.to_vtk(fname='beforeopt_EPI.vtu')
 mesh.delete_degree_4()
 isBdNode = mesh.ds.boundary_node_flag()
 isFreeNode = ~isBdNode
+
+node = mesh.entity('node')
+cell = mesh.entity('cell')
 # 去除度为4的点后进行优化
-opt = TriRadiusRatio(mesh)
-opt.iterate_solver(method='Bjacobi',isFreeNode=isFreeNode)
+#opt = TriRadiusRatio(mesh)
+
+from app.HVSCMesh.optimizer import *
+from fealpy.experimental.mesh import TriangleMesh
+
+mesh = TriangleMesh(node,cell)
+mesh_quality = RadiusRatioQuality(mesh)
+node = mesh.entity('node')
+mesh = iterate_solver(mesh)
+
 angle = mesh.angle()
 max_angle = np.max(angle,axis=1)
 angles = max_angle*(180/np.pi)
@@ -38,10 +50,9 @@
 fig,axes1= plt.subplots()
 axes1.set_ylim(0,8500)
 fig.text(0.15,0.85,f'Number of angles>90:{m90}',fontsize=12)
-mesh.show_angle(axes1,max_angle)
+show_angle(axes1,max_angle)
 plt.show()
 
 mesh.celldata['angles'] = angles
 mesh.to_vtk(fname='optEPI.vtu')
-gmsh.finalize()
 

app/HVSCMesh/area_adaptive/rreo_RB.py

@@ -7,20 +7,6 @@
 
 from fealpy.mesh import TriangleMesh
 
-'''
-gmsh.initialize()
-gmsh.merge('adaptive_case_4_RB_IGCT.msh')
-
-ntags, vxyz, _ = gmsh.model.mesh.getNodes()
-node = vxyz.reshape((-1,3))
-node = node[:,:2]
-vmap = dict({j:i for i,j in enumerate(ntags)})
-tris_tags,evtags = gmsh.model.mesh.getElementsByType(2)
-evid = np.array([vmap[j] for j in evtags])
-cell = evid.reshape((tris_tags.shape[-1],-1))
-mesh = TriangleMesh(node,cell)
-'''
-
 mesh = TriangleMesh.from_meshio('area_RB_adaptive.vtu')
 node = mesh.entity('node')
 cell = mesh.entity('cell')
@@ -48,7 +34,9 @@
 # 去除度为4的点后进行优化
 #opt = TriRadiusRatio(mesh)
 
-from optimizer import *
+from app.HVSCMesh.optimizer import *
+from fealpy.experimental.mesh import TriangleMesh
+
 mesh = TriangleMesh(node,cell)
 mesh_quality = RadiusRatioQuality(mesh)
 node = mesh.entity('node')

app/HVSCMesh/optimizer.py

@@ -6,7 +6,7 @@
 from fealpy.experimental.mesh.triangle_mesh import TriangleMesh
 from fealpy.experimental.mesh.tetrahedron_mesh import TetrahedronMesh
 from fealpy.experimental.mesh.mesh_quality import RadiusRatioQuality
-from radius_ratio_objective import RadiusRatioSumObjective
+from .radius_ratio_objective import RadiusRatioSumObjective
 
 def show_mesh_quality(q1,ylim=1000):
     fig,axes= plt.subplots()
@@ -52,6 +52,31 @@ def show_mesh_quality(q1,ylim=1000):
     plt.show()
     return 0
 
+def show_angle(axes, angle):
+        """
+        @brief 显示网格角度的分布直方图
+        """
+        hist, bins = np.histogram(angle.flatten('F') * 180 / bm.pi, bins=50, range=(0, 180))
+        center = (bins[:-1] + bins[1:]) / 2
+        axes.bar(center, hist, align='center', width=180 / 50.0)
+        axes.set_xlim(0, 180)
+        mina = np.min(angle.flatten('F') * 180 / np.pi)
+        maxa = np.max(angle.flatten('F') * 180 / np.pi)
+        meana = np.mean(angle.flatten('F') * 180 / np.pi)
+        axes.annotate('Min angle: {:.4}'.format(mina), xy=(0.41, 0.5),
+                      xytext=(0.15,0.85),
+                      textcoords='axes fraction',
+                      horizontalalignment='left', verticalalignment='top')
+        axes.annotate('Max angle: {:.4}'.format(maxa), xy=(0.41, 0.45),
+                      xytext=(0.15,0.75),
+                      textcoords='axes fraction',
+                      horizontalalignment='left', verticalalignment='top')
+        axes.annotate('Average angle: {:.4}'.format(meana), xy=(0.41, 0.40),
+                      xytext=(0.15,0.7),
+                      textcoords='axes fraction',
+                      horizontalalignment='left', verticalalignment='top')
+        return mina, maxa, meana
+
 def BlockJacobi2d(node,A,B,isFreeNode):
     NN = node.shape[0] 
     isBdNode = ~isFreeNode

fealpy/experimental/mesh/triangle_mesh.py

@@ -335,7 +335,16 @@ def circumcenter(self, index: Index=_S, returnradius=False):
             return c
 
     def angle(self):
-        pass
+        NC = self.number_of_cells()
+        cell = self.entity('cell')
+        node = self.entity('node')
+        localEdge = self.localEdge
+        angle = bm.zeros((NC, 3), dtype=self.ftype)
+        for i,(j,k) in zip(range(3),localEdge):
+            v0 = node[cell[:, j]] - node[cell[:, i]]
+            v1 = node[cell[:, k]] - node[cell[:, i]]
+            angle[:,i] = bm.arccos(bm.sum(v0*v1,axis=1)/bm.sqrt(bm.sum(v0**2,axis=1)*bm.sum(v1**2,axis=1)))
+        return angle
 
     def show_angle(self, axes, angle=None):
         """