From 5fb4ff66deb0ef4f5b779795cb2447b285c4149d Mon Sep 17 00:00:00 2001 From: Christophe Bourcier Date: Fri, 21 Apr 2023 14:27:50 +0200 Subject: [PATCH] Add examples to salome test. And add some checks to test the generated mesh --- doc/salome/examples/examples.set | 6 +- doc/salome/examples/ghs3dSetParametersDemo.py | 8 +++ doc/salome/examples/ghs3d_enfmesh.py | 57 +++++++++++++++++-- doc/salome/examples/ghs3d_enfvert.py | 47 ++++++++++++--- doc/salome/examples/ghs3d_optimization.py | 20 +++++-- doc/salome/examples/ghs3ddemo.py | 10 +++- 6 files changed, 129 insertions(+), 19 deletions(-) diff --git a/doc/salome/examples/examples.set b/doc/salome/examples/examples.set index 11f137d..ce85751 100644 --- a/doc/salome/examples/examples.set +++ b/doc/salome/examples/examples.set @@ -19,4 +19,8 @@ SET(EXAMPLE_NAMES ghs3dSetParametersDemo -) \ No newline at end of file + ghs3ddemo + ghs3d_enfmesh + ghs3d_enfvert + ghs3d_optimization +) diff --git a/doc/salome/examples/ghs3dSetParametersDemo.py b/doc/salome/examples/ghs3dSetParametersDemo.py index b2c6f4c..82afb52 100644 --- a/doc/salome/examples/ghs3dSetParametersDemo.py +++ b/doc/salome/examples/ghs3dSetParametersDemo.py @@ -13,6 +13,8 @@ smesh = smeshBuilder.New() box = geompy.MakeBoxDXDYDZ(200., 200., 200.) geompy.addToStudy(box, "box") +expected_volume = 200**3 + # create a mesh on the box mgtetraMesh = smesh.Mesh(box,"box: MG-Tetra and NETGEN_1D_2D mesh") @@ -30,6 +32,9 @@ MG_Tetra_Parameters_1.SetPthreadMode( 1 ) # 0 - none, 1 - aggressive, 2 - status = mgtetraMesh.Compute() assert( status ) +volume = smesh.GetVolume(mgtetraMesh) +assert (volume-expected_volume)/expected_volume < 1e-12 + mgtetraHPCMesh = smesh.Mesh(box,"box: MG-Tetra HPC and NETGEN_1D_2D mesh") status = mgtetraHPCMesh.AddHypothesis(NETGEN_2D_Parameters_1) NETGEN_1D_2D_1 = mgtetraHPCMesh.Triangle(algo=smeshBuilder.NETGEN_1D2D) @@ -44,6 +49,9 @@ MG_Tetra_Parameters_2.SetParallelMode( 1 ) # 0 - none, 1 - reproducible_giv status = mgtetraHPCMesh.Compute() assert( status ) +volume_HPC = smesh.GetVolume(mgtetraHPCMesh) +assert (volume_HPC-expected_volume)/expected_volume < 1e-12 + # End of script diff --git a/doc/salome/examples/ghs3d_enfmesh.py b/doc/salome/examples/ghs3d_enfmesh.py index 9752289..d05d718 100644 --- a/doc/salome/examples/ghs3d_enfmesh.py +++ b/doc/salome/examples/ghs3d_enfmesh.py @@ -62,13 +62,18 @@ Mesh_cylindre.AddHypothesis( MG_CADSurf_Parameters2 ) face_cyl_faces = Mesh_cylindre.GroupOnGeom(face_cyl,'group_face_cyl', SMESH.FACE) face_cyl_edges = Mesh_cylindre.GroupOnGeom(face_cyl,'group_edge_cyl', SMESH.EDGE) face_cyl_nodes = Mesh_cylindre.GroupOnGeom(face_cyl,'group_node_cyl', SMESH.NODE) -Mesh_cylindre.Compute() +ok = Mesh_cylindre.Compute() +if not ok: + raise Exception("Error when computing Mesh_cylindre with MG_CADSurf") # Create the mesh on the cylinder Mesh_box_tri = smesh.Mesh(box,"Mesh_box_tri") Mesh_box_tri.AddHypothesis( MG_CADSurf ) Mesh_box_tri.AddHypothesis( MG_CADSurf_Parameters ) -Mesh_box_tri.Compute() +ok = Mesh_box_tri.Compute() +if not ok: + raise Exception("Error when computing Mesh_box_tri with MG_CADSurf") + # Create 4 copies of the 2D mesh to test the 3 types of contraints (NODE, EDGE, FACE) # from the whole mesh and from groups of elements. @@ -97,9 +102,49 @@ MG_Tetra_Parameters_edge.SetEnforcedMeshWithGroup(face_cyl_edges,SMESH.EDGE,"edg MG_Tetra_Parameters_face.SetEnforcedMeshWithGroup(face_cyl_faces,SMESH.FACE,"faces from face_cyl_faces") #Compute the meshes -mesh_node.Compute() -mesh_edge.Compute() -mesh_face.Compute() -mesh_mesh.Compute() +ok = mesh_node.Compute() +if not ok: + raise Exception("Error when computing mesh_node") +ok = mesh_edge.Compute() +if not ok: + raise Exception("Error when computing mesh_edge") +ok = mesh_face.Compute() +if not ok: + raise Exception("Error when computing mesh_face") +ok = mesh_mesh.Compute() +if not ok: + raise Exception("Error when computing mesh_mesh") + +## Get the number of nodes of a mesh or a group of mesh +def getNbNodes(mesh): + try: + #mesh + nb_nodes = mesh.NbNodes() + except: + # group + nb_nodes = mesh.GetNumberOfNodes() + return nb_nodes + +## Check enforced mesh is correct +def checkEnforcedMesh(mesh, enforced_group, result_group, tol=1e-7): + name = "Check " + mesh.GetName() + # test nb nodes are the same between enforced group and the result group + if getNbNodes(enforced_group) != getNbNodes(enforced_group): + raise Exception("Enforce mesh failed for %s, wrong nummber of nodes"%name) + # test that all nodes are imposed at the same place by finding coincident nodes in a compound mesh + mesh_check = smesh.Concatenate([enforced_group, result_group], 0, name=name) + ll_coincident_nodes = mesh_check.FindCoincidentNodes(tol) + coincident_nodes = [item for sublist in ll_coincident_nodes for item in sublist] + mesh_check.MakeGroupByIds("coincident_nodes", SMESH.NODE, coincident_nodes) + mesh_nodes = mesh_check.GetNodesId() + if len(ll_coincident_nodes) != getNbNodes(enforced_group): + non_coincident_nodes = list(set(mesh_nodes) - set(coincident_nodes)) + mesh_check.MakeGroupByIds("non_coincident_nodes", SMESH.NODE, non_coincident_nodes) + raise Exception("Enforce mesh failed for %s"%name) + +checkEnforcedMesh(mesh_node, face_cyl_nodes, mesh_node.GetGroupByName("nodes from face_cyl_nodes")[0]) +checkEnforcedMesh(mesh_edge, face_cyl_edges, mesh_edge.GetGroupByName("edges from face_cyl_edges")[0]) +checkEnforcedMesh(mesh_face, face_cyl_faces, mesh_face.GetGroupByName("faces from face_cyl_faces")[0]) +checkEnforcedMesh(mesh_mesh, Mesh_cylindre.GetMesh(), mesh_mesh.GetGroupByName("faces from cylinder")[0]) # End of script diff --git a/doc/salome/examples/ghs3d_enfvert.py b/doc/salome/examples/ghs3d_enfvert.py index 14c25cd..c4cf165 100644 --- a/doc/salome/examples/ghs3d_enfvert.py +++ b/doc/salome/examples/ghs3d_enfvert.py @@ -10,6 +10,8 @@ # Ex1: Add one enforced vertex with coordinates (50,50,100) # and physical size 2. +import math + import salome salome.salome_init() @@ -20,6 +22,13 @@ import SMESH from salome.smesh import smeshBuilder smesh = smeshBuilder.New() +## check if a node is at the coords given the tolerance +def isNodeAtCoord(node_id, x, y, z, tol=1e-12): + xn, yn, zn = mgtetraMesh.GetNodeXYZ(node_id) + dist = math.sqrt((x-xn)**2+(y-yn)**2+(z-zn)**2) + print("dist: ", dist) + return dist < tol + # create a box box = geompy.MakeBoxDXDYDZ(200., 200., 200.) geompy.addToStudy(box, "box") @@ -28,7 +37,9 @@ mgtetraMesh = smesh.Mesh(box,"box: MG-Tetra and MG-CADSurf mesh") # create a MG-CADSurf algorithm for faces mgtetraMesh.Triangle(algo=smeshBuilder.MG_CADSurf) # compute the mesh -mgtetraMesh.Compute() +ok = mgtetraMesh.Compute() +if not ok: + raise Exception("Error when computing MG-CADSurf mesh") # Make a copy of the 2D mesh mgtetraMesh_wo_geometry = smesh.CopyMesh( mgtetraMesh, 'MG-Tetra w/o geometry', 0, 0) @@ -37,10 +48,19 @@ mgtetraMesh_wo_geometry = smesh.CopyMesh( mgtetraMesh, 'MG-Tetra w/o geometry', MG_Tetra = mgtetraMesh.Tetrahedron( smeshBuilder.MG_Tetra ) MG_Tetra_Parameters = MG_Tetra.Parameters() # Create the enforced vertex -MG_Tetra_Parameters.SetEnforcedVertex( 50, 50, 100, 2) # no group +x1 = 50 +y1 = 50 +z1 = 100 +MG_Tetra_Parameters.SetEnforcedVertex( x1, y1, z1, 2) # no group # Compute the mesh -mgtetraMesh.Compute() +ok = mgtetraMesh.Compute() +if not ok: + raise Exception("Error when computing MG_Tetra mesh") + +# Check that the enforced node is at the enforced coords +node_closest = mgtetraMesh.FindNodeClosestTo(x1, y1, z1) +assert isNodeAtCoord(node_closest, x1, y1, z1) # Ex2: Add one vertex enforced by a GEOM vertex at (50,50,100) # with physical size 5 and add it to a group called "My special nodes" @@ -51,13 +71,26 @@ mgtetraMesh_wo_geometry.AddHypothesis( MG_Tetra ) mgtetraMesh_wo_geometry.AddHypothesis( MG_Tetra_Parameters_wo_geometry ) # Create the enforced vertex -p1 = geompy.MakeVertex(150, 150, 100) -geompy.addToStudy(p1, "p1") -MG_Tetra_Parameters_wo_geometry.SetEnforcedVertexGeomWithGroup( p1, 5 , "My special nodes") +x2 = 50 +y2 = 50 +z2 = 100 +p2 = geompy.MakeVertex(x2, y2, z2) +geompy.addToStudy(p2, "p2") +gr_enforced_name = "My special nodes" +MG_Tetra_Parameters_wo_geometry.SetEnforcedVertexGeomWithGroup( p2, 5 , gr_enforced_name) #MG_Tetra_Parameters.SetEnforcedVertexGeom( p1, 5 ) # no group # compute the mesh -mgtetraMesh_wo_geometry.Compute() +ok = mgtetraMesh_wo_geometry.Compute() +if not ok: + raise Exception("Error when computing MG_Tetra mesh without geometry") + +# Check that the enforced node is at the enforced coords +gr_enforced_nodes = mgtetraMesh_wo_geometry.GetGroupByName(gr_enforced_name)[0] +assert (gr_enforced_nodes.Size() == 1) +node_enforced = gr_enforced_nodes.GetIDs()[0] + +assert isNodeAtCoord(node_enforced, x2, y2, z2) # Erase all enforced vertices MG_Tetra_Parameters.ClearEnforcedVertices() diff --git a/doc/salome/examples/ghs3d_optimization.py b/doc/salome/examples/ghs3d_optimization.py index 34d9e6d..baa6d7f 100644 --- a/doc/salome/examples/ghs3d_optimization.py +++ b/doc/salome/examples/ghs3d_optimization.py @@ -2,10 +2,11 @@ import salome salome.salome_init() from salome.geom import geomBuilder -geompy = geomBuilder.New(salome.myStudy) +geompy = geomBuilder.New() +import SMESH from salome.smesh import smeshBuilder -smesh = smeshBuilder.New(salome.myStudy) +smesh = smeshBuilder.New() # create a disk disk = geompy.MakeDiskR(100., 1, theName="disk") @@ -14,7 +15,9 @@ disk = geompy.MakeDiskR(100., 1, theName="disk") mesh = smesh.Mesh( disk ) cadsurf = mesh.Triangle( smeshBuilder.MG_CADSurf ) cadsurf.SetQuadAllowed( True ) -mesh.Compute() +ok = mesh.Compute() +if not ok: + raise Exception("Error when computing mesh") # extrude the 2D mesh into a prismatic mesh mesh.ExtrusionSweepObject( mesh, [0,0,10], 7 ) @@ -32,10 +35,19 @@ mg_opt.SetSmoothOffSlivers( True ) mg_opt.SetOptimizationLevel( smeshBuilder.Strong_Optimization ) # run optimization -optMesh.Compute() +ok = optMesh.Compute() +if not ok: + raise Exception("Error when computing optimization mesh") print("Nb tetra before optimization", mesh.NbTetras()) print("Nb tetra after optimization", optMesh.NbTetras()) +# Check that aspect ratio 3D of optimized mesh is better than original mesh +min_aspectratio_orig, max_aspectratio_orig = mesh.GetMinMax(SMESH.FT_AspectRatio3D) +min_aspectratio_optim, max_aspectratio_optim = optMesh.GetMinMax(SMESH.FT_AspectRatio3D) + +assert (min_aspectratio_orig - min_aspectratio_optim)/min_aspectratio_orig < 0.5 +assert (max_aspectratio_orig - max_aspectratio_optim)/max_aspectratio_orig < 0.5 + # End of script diff --git a/doc/salome/examples/ghs3ddemo.py b/doc/salome/examples/ghs3ddemo.py index 46d098a..c258692 100644 --- a/doc/salome/examples/ghs3ddemo.py +++ b/doc/salome/examples/ghs3ddemo.py @@ -19,7 +19,15 @@ MG_CADSurf = mgtetraMesh.Triangle(algo=smeshBuilder.MG_CADSurf) MG_Tetra = mgtetraMesh.Tetrahedron(algo=smeshBuilder.MG_Tetra) # compute the mesh -mgtetraMesh.Compute() +ok = mgtetraMesh.Compute() + +if not ok: + raise Exception("Error when computing mgtetraMesh") + +volume = smesh.GetVolume(mgtetraMesh) +expected_volume = 200**3 + +assert (volume-expected_volume)/expected_volume < 1e-12 # End of script -- 2.39.2