doc/salome/examples/curvature_face.py

   1 # Curvature of a Face along given direction
   2
   3 import salome
   4 salome.salome_init_without_session()
   5 import GEOM
   6 from salome.geom import geomBuilder
   7 geompy = geomBuilder.New()
   8 import math
   9 import numpy as np
  10
  11 def test_acceptance():
  12   """
  13   Acceptance test [tuleap29472]
  14   """
  15   Vector = [0,100,100]
  16   O = geompy.MakeVertex(0, 0, 0)
  17   OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
  18   OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
  19   OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
  20   Cylinder_1 = geompy.MakeCylinderRH(100, 300)
  21   Translation_1 = geompy.MakeTranslation(Cylinder_1, 0, 0, -150)
  22   Vertex_1 = geompy.MakeVertex(100, 0, 0)
  23   Vertex_2 = geompy.MakeVertex(100, -Vector[2], Vector[1])
  24   Line_1 = geompy.MakeLineTwoPnt(Vertex_1, Vertex_2)
  25   Plane_1 = geompy.MakePlane(Vertex_1, Line_1, 2000)
  26   Rotation_1 = geompy.MakeRotation(Translation_1, OZ, 90*math.pi/180.0)# avoid to have degenerated edge across Vertex_1
  27
  28   [Face_1,Face_2,Face_3] = geompy.ExtractShapes(Rotation_1, geompy.ShapeType["FACE"], True)
  29
  30   curvature_29472 = np.array( geompy.VectorCoordinates( geompy.CurvatureOnFace(Face_2, Vertex_1, geompy.MakeVectorDXDYDZ(*Vector))) ).reshape(1,3)
  31   expected_curvature = np.array( [-200.0,0.0,0.0] ).reshape(1,3)
  32   assert( np.isclose( 0.0, np.linalg.norm( curvature_29472 - expected_curvature ) ,rtol=0,atol=1e-5 ) )
  33
  34   Intersection_1 = geompy.MakeSection(Face_2, Plane_1, True)
  35   geompy.addToStudy( O, 'O' )
  36   geompy.addToStudy( OX, 'OX' )
  37   geompy.addToStudy( OY, 'OY' )
  38   geompy.addToStudy( OZ, 'OZ' )
  39   geompy.addToStudy( Vertex_1, 'Vertex_1' )
  40   geompy.addToStudy( Cylinder_1, 'Cylinder_1' )
  41   geompy.addToStudy( Translation_1, 'Translation_1' )
  42   geompy.addToStudy( Vertex_2, 'Vertex_2' )
  43   geompy.addToStudy( Line_1, 'Line_1' )
  44   geompy.addToStudy( Plane_1, 'Plane_1' )
  45   geompy.addToStudy( Rotation_1, 'Rotation_1' )
  46   geompy.addToStudyInFather( Rotation_1, Face_1, 'Face_1' )
  47   geompy.addToStudyInFather( Rotation_1, Face_2, 'Face_2' )
  48   geompy.addToStudyInFather( Rotation_1, Face_3, 'Face_3' )
  49   geompy.addToStudy( Intersection_1, 'Intersection_1' )
  50   angle = math.asin(Vector[2]/math.sqrt(Vector[1]*Vector[1]+Vector[2]*Vector[2]))
  51   tmp = geompy.MakeTranslation(Intersection_1,*[-elt for elt in geompy.PointCoordinates(Vertex_1)])
  52   tmp = geompy.MakeRotation(tmp,OX,-angle)
  53   Intersection_1_OXY = geompy.MakeTranslation(tmp,*geompy.PointCoordinates(Vertex_1))
  54   geompy.addToStudy( Intersection_1_OXY, 'Intersection_1_OXY' )
  55
  56   eps = 0.01
  57   offset = 0.75
  58   p0 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset-eps) ) ).reshape(1,3)
  59   p1 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset) ) ).reshape(1,3)
  60   p2 = np.array( geompy.PointCoordinates( geompy.MakeVertexOnCurve(Intersection_1_OXY,offset+eps) ) ).reshape(1,3)
  61   assert( np.isclose(0.0,np.linalg.norm(p1- np.array(geompy.PointCoordinates(Vertex_1)).reshape(1,3)  ),rtol=0,atol=1e-8) )
  62   p01=(p0+p1)/2
  63   p12=(p1+p2)/2
  64   v0 = (p1-p0)/np.linalg.norm(p1-p0)
  65   v1 = (p2-p1)/np.linalg.norm(p2-p1)
  66   computedRadius =  1/np.linalg.norm((v1-v0)/np.linalg.norm(p12-p01))
  67   # manual detection of radius : https://fr.wikipedia.org/wiki/Courbure_d%27un_arc
  68   circle = geompy.MakeCircle(O,OZ,computedRadius)
  69   circle = geompy.MakeTranslation(circle,100-computedRadius,0,0)
  70   geompy.addToStudy(circle, "expectedCircle")
  71   print("Radius expected is {}".format(computedRadius))
  72   print("Radius obtain by CurvatureOnFace is {}".format(np.linalg.norm(curvature_29472)))
  73
  74 O = geompy.MakeVertex(0, 0, 0, 'O')
  75 OX = geompy.MakeVectorDXDYDZ(1, 0, 0, 'OX')
  76 OY = geompy.MakeVectorDXDYDZ(0, 1, 0, 'OY')
  77 OZ = geompy.MakeVectorDXDYDZ(0, 0, 1, 'OZ')
  78
  79 pXYZ = geompy.MakeVertex(105, 105, 105, 'pXYZ')
  80 pY = geompy.MakeVertex(0, 105, 0, 'pY')
  81 pZ = geompy.MakeVertex(0, 0, 105, 'pZ')
  82
  83 vZ_XY = geompy.MakeVectorDXDYDZ(-1, -1, 1, 'vZ-XY')
  84 vZ_XY2 = geompy.MakeVectorDXDYDZ(-1, -1, 10, 'vZ-XY')
  85 vZ_XY3 = geompy.MakeVectorDXDYDZ(-1, -1, 100, 'vZ-XY')
  86
  87 R = 100.0
  88
  89 # I. Curvature of a Sphere
  90 Sphere_1 = geompy.MakeSphereR(R, 'Sphere_1')
  91 [Sph] = geompy.ExtractShapes(Sphere_1, geompy.ShapeType["FACE"], True, "Sph")
  92
  93 curvature_1 = geompy.CurvatureOnFace(Sph, pXYZ, OX,    'curvature_sph_pXYZ_OX')
  94 curvature_2 = geompy.CurvatureOnFace(Sph, pXYZ, vZ_XY, 'curvature_sph_pXYZ_vt')
  95 curvature_3 = geompy.CurvatureOnFace(Sph, pY,   OX,    'curvature_sph_pY_OX')
  96
  97 # All sphere curvature radiuces = R
  98 assert(abs(geompy.BasicProperties(curvature_1)[0] - R) < 1e-07)
  99 assert(abs(geompy.BasicProperties(curvature_2)[0] - R) < 1e-07)
 100 assert(abs(geompy.BasicProperties(curvature_3)[0] - R) < 1e-07)
 101
 102 # Pole
 103 isExcept = False
 104 try:
 105   geompy.CurvatureOnFace(Sph, pZ, OX)
 106 except:
 107   isExcept = True
 108 assert(isExcept)
 109
 110 # Normal direction
 111 isExcept = False
 112 try:
 113   geompy.CurvatureOnFace(Sph, pY,   OY)
 114 except:
 115   isExcept = True
 116 assert(isExcept)
 117
 118 # II. Curvature of a Cylinder
 119 Cylinder_1 = geompy.MakeCylinderRH(R, 300, 'Cylinder_1')
 120 [Face_1,Face_2,Face_3] = geompy.ExtractShapes(Cylinder_1, geompy.ShapeType["FACE"], True, "Face")
 121
 122 # Curvature radius of a cylinder along any direction, orthogonal to its Z axis, equal to R
 123 curvature_4 = geompy.CurvatureOnFace(Face_2, pY, OX, 'curvature_cyl_pY_OX')
 124 assert(abs(geompy.BasicProperties(curvature_4)[0] - R) < 1e-07)
 125
 126 # Curvature radius of a cylinder along its Z direction is infinite
 127 curvature_zero = geompy.CurvatureOnFace(Face_2, pY, OZ)
 128 assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
 129 assert(not curvature_zero)
 130
 131 # Curvature radius of a cylinder along some direction, different from two above
 132 curvature_5 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY,  'curvature_cyl_pY_vZ_XY')
 133 curvature_6 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY2, 'curvature_cyl_pY_vZ_XY2')
 134 curvature_7 = geompy.CurvatureOnFace(Face_2, pY, vZ_XY3, 'curvature_cyl_pY_vZ_XY3')
 135
 136 # R < r5 < r6 < r7
 137 # r5 = 100.01, r6 = 101.0, r7 = 200
 138 r5 = geompy.BasicProperties(curvature_5)[0]
 139 r6 = geompy.BasicProperties(curvature_6)[0]
 140 r7 = geompy.BasicProperties(curvature_7)[0]
 141
 142 assert(R  + 1e-07 < r5)
 143 assert(r5 + 1e-07 < r6)
 144 assert(r6 + 1e-07 < r7)
 145
 146 # Projection aborted. Point is out of the face boundaries.
 147 isExcept = False
 148 try:
 149   pXY_Z = geompy.MakeVertex(105, 105, -105, 'pXY_Z')
 150   geompy.CurvatureOnFace(Face_2, pXY_Z, OX, 'curvature_cyl_pXY_Z')
 151 except:
 152   isExcept = True
 153 assert(isExcept)
 154
 155 # Projection aborted (point on axis). Equal distances to many points.
 156 isExcept = False
 157 try:
 158   geompy.CurvatureOnFace(Face_2, O, vZ_XY, 'curvature_cyl_O')
 159 except:
 160   isExcept = True
 161 assert(isExcept)
 162
 163 # Curvature radius of a planar face is infinite
 164 curvature_zero_2 = geompy.CurvatureOnFace(Face_1, pZ, OX)
 165 assert(geompy.MeasuOp.GetErrorCode() == "ZERO_CURVATURE")
 166 assert(not curvature_zero_2)
 167
 168 # III. Curvature of a "Horse saddle"
 169 [Edge_1,Edge_2,Edge_3] = geompy.ExtractShapes(Sphere_1, geompy.ShapeType["EDGE"], True)
 170 geompy.addToStudyInFather( Sphere_1, Edge_1, 'Edge_1' )
 171 geompy.addToStudyInFather( Sphere_1, Edge_2, 'Edge_2' )
 172 geompy.addToStudyInFather( Sphere_1, Edge_3, 'Edge_3' )
 173
 174 Rotation_1 = geompy.MakeRotation(Edge_3, OX, 90*math.pi/180.0, 'Rotation_1')
 175 Rotation_2 = geompy.MakeRotation(Rotation_1, OY, 180*math.pi/180.0, 'Rotation_2')
 176 Translation_1 = geompy.MakeTranslation(Rotation_2, 200, 0, 0, 'Translation_1')
 177 Translation_2 = geompy.MakeTranslation(Edge_3, 100, 100, 0, 'Translation_2')
 178 Translation_3 = geompy.MakeTranslation(Translation_2, 0, -200, 0, 'Translation_3')
 179 Filling_1 = geompy.MakeFilling([Translation_2, Edge_3, Translation_3])
 180 geompy.addToStudy(Filling_1, 'Filling_1')
 181 Vertex_2 = geompy.MakeVertex(100, 0, 0, 'Vertex_2')
 182
 183 curvature_Y = geompy.CurvatureOnFace(Filling_1, Vertex_2, OY, 'curvature_Y')
 184 curvature_Z = geompy.CurvatureOnFace(Filling_1, Vertex_2, OZ, 'curvature_Z')
 185
 186 cury = np.array( geompy.VectorCoordinates(curvature_Y) ).reshape(1,3)
 187 curz = np.array( geompy.VectorCoordinates(curvature_Z) ).reshape(1,3)
 188 cury_expected = np.array( [50,0,0] ).reshape(1,3)
 189 curz_expected = np.array( [-100,0,0] ).reshape(1,3)
 190 assert( np.isclose( 0.0, np.linalg.norm( cury - cury_expected ) ,rtol=0,atol=1e-5 ) )
 191 assert( np.isclose( 0.0, np.linalg.norm( curz - curz_expected ) ,rtol=0,atol=1e-5 ) )
 192
 193 # Normal direction
 194 norm_1 = geompy.GetNormal(Filling_1, Vertex_2, "Normal_1")
 195 isExcept = False
 196 try:
 197   geompy.CurvatureOnFace(Filling_1, Vertex_2, norm_1)
 198 except:
 199   isExcept = True
 200 assert(isExcept)
 201
 202 # acceptance case
 203 test_acceptance()