ochubar / Radia

3D Magnetostatics Computer Code
Other
35 stars 20 forks source link

ObjPolyhdr - wrong coordinate #16

Open quorx opened 3 years ago

quorx commented 3 years ago

Hi,

Possible issue with the ObjPolyhdr function. Here's a simple example that demonstrates the issue for a single tetrahedral element.

import radia as rad 

rad.UtiDelAll()

zPnt = 0.249819

p1 = [1.00000,       0.26250,       0.25000]
p2 = [0.80000,       0.00000,       0.25000]
p3 = [0.80000,       0.25000,       0.50000]
p4 = [0.50000,       0.27500,          zPnt]

connectivity = [[1,3,2],[1,2,4],[2,3,4],[3,1,4]]

tet1 = rad.ObjPolyhdr([p1, p2, p3, p4], connectivity)

print(rad.UtiDmp(tet1))

Which gives the following listing of the element

Index 1: Magnetic field source object: Relaxable: Polyhedron
   Number of faces: 4
   {x,y,z}= {0.775,0.196875,0.312455}
   {mx,my,mz}= {0,0,0}
   Face Vertices:
   {{1,0.2625,0.25},{0.8,0.25,0.5},{0.8,0,0.25}},
   {{1,0.2625,0.249819},{0.8,0,0.249819},{0.5,0.275,0.249819}},
   {{0.8,2.77556e-17,0.25},{0.8,0.25,0.5},{0.5,0.275,0.249819}},
   {{0.8,0.25,0.5},{1,0.2625,0.25},{0.5,0.275,0.249819}},

   Material applied: None
   Transformations applied: None
   Memory occupied: 1936 bytes

The second face connects points p1, p2 & p4, but the z-coordinate of p4 has overwritten the z-coordinates of p1 and p2. If the value zPnt is made smaller (0.2498 for example), the problem goes away. So it's only an error of order 1 part in 1000, which remains true if all coordinates are scaled.

The error is small so probably won't affect field results much, but when importing a tet mesh into Radia and then exporting the model and results (via ObjDrwVTK) for display in Paraview (after culling interior faces) there are interior faces visible that shouldn't be (see attached images with and without internal facets).

Thank you!

tet mesh interior facets tet mesh without interior facets