Current ouput .out reports nor SurfaceCurrentsAndCharges Dataset not matching POSTFEKO results

Ricardo Sendrea
Ricardo Sendrea Altair Community Member
edited August 8 in Community Q&A

Hello,

There may be a bug in how the currents are stored (or ordered) in the .out files. Using the Surface Current and Charges dataset object within the LUA environment, we wrote the results into a text file (please see the script below for details). We tested this with purely dielectric (epsilon_r = 2.56) and purely metallic sphere models (radius = 0.5m). Here, the spheres are illuminated by a plane wave and solved at 380MHz. As you can see from the figure below (FEKO plots on the left and the MATLAB via textfile export on the right) the current magnitudes are slightly off from FEKO, but the current distributions are correct.

image

We repeated this with a coated sphere model (radius = 0.5m), where the interior is a metal sphere (radius = 0.25m) enveloped by a dielectric sphere (epsilon_r = 2.56). Here, it is illuminated by a plane wave solved at 250MHz. Unlike the other cases, we noticed that even though the electric currents along the interior sphere are identical to FEKO, the electric and magnetic currents on the dielectric region are completely wrong. As shown in the figure below, the electric currents along the metallic region are as expected from FEKO, but the electric currents along the dielectric (although nearly the expected magnitude) are completely wrong. We cross-validated our textfile with the outputs from the CADFEKO .out and  POSTFEKO .os / .ol files and they are the same.

image

I suspect there may be an ordering (i.e., indexing) problem, for example, the current values are stored at the wrong triangle index. Regarding the difference in magnitude for the purely conductive and dielectric sphere cases, this doesn't occur if you output the currents via the LUA DataSet: ExportMatFile() function, which unfortunately doesn't work for composite designs in FEKO 2020.0. Perhaps this is a consequence of the same bug?

The LUA code:

app = pf.GetApplication() app:NewProject() app:OpenFile([[path_here\coated_sphere.fek]]) -- Set up the file details absPath  = [[path_here\]]  -- HIGHLY recommended to hard code the absolute path! meshfile = "coated_sphere_sol.txt" -- Find all possible data files names = pf.SurfaceCurrentsAndCharges.GetNames() --Get characteristic mode model currents currents = pf.SurfaceCurrentsAndCharges.GetDataSet(names[1]) --Load a given model sConf = app.Models["coated_sphere"].Configurations[1] --Import the mesh mesh = sConf.Mesh count = mesh.TriangleFaces.Count --Number of Triangles --Get mesh values points = mesh.Points --I/O for mesh write upper-- Open a file to write to. file = io.open(absPath..meshfile, "w") --Create Column Headers Vertice Index - Position(X,Y,Z) colHeader = string.format("%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s\n", "Frequency", "Element","X", "Y", "Z", "real_Jx", "realJy", "real_Jz", "imag_Jx", "imag_Jy", "imag_Jz", "real_Mx", "real_My", "real_Mz", "imag_Mx", "imag_My", "imag_Mz") --Write Column header to file file:write(colHeader) -- Get vector of all solved frequencies frequencies = currents.Axes[1] for c = 1, #frequencies do         str = ""         elementNum = 1         pointCount = 1         for i = 1, count do             --Load mesh triangle             triangle = mesh.TriangleFaces[i].Triangles             count2 = triangle.Count --Number of triangle vertices               for j = 1, count2 do             info = triangle[j].VertexIndices                                          -- Pre-Allocate Center of Triangle Value for ouput              -- Electric Current             real_Js_x = 0             real_Js_y = 0             real_Js_z = 0             imag_Js_x = 0             imag_Js_y = 0             imag_Js_z = 0                          -- Magnetic Current             real_Ms_x = 0             real_Ms_y = 0             real_Ms_z = 0             imag_Ms_x = 0             imag_Ms_y = 0             imag_Ms_z = 0                          -- Element Centroid             centroidx = 0             centroidy = 0             centroidz = 0                          for ii = 1, #info do                                                    -- Get vertice location                 position = points[triangle[j].VertexIndices[ii]]                 -- Add to centroid                 centroidx = centroidx + position[1]                 centroidy = centroidy + position[2]                 centroidz = centroidz + position[3]                                  -- Add Electric Currents                 real_Js_x = currents[c][pointCount].ElectricX:real()                 real_Js_y = currents[c][pointCount].ElectricY:real()                 real_Js_z = currents[c][pointCount].ElectricZ:real()                                 -- Imag component                 imag_Js_x = currents[c][pointCount].ElectricX:imag()                 imag_Js_y = currents[c][pointCount].ElectricY:imag()                 imag_Js_z = currents[c][pointCount].ElectricZ:imag()                                      if currents[c][pointCount].MagneticX ~= nil then                 -- Add Magnetic Currents                 real_Ms_x = currents[c][pointCount].MagneticX:real()                 real_Ms_y = currents[c][pointCount].MagneticY:real()                 real_Ms_z = currents[c][pointCount].MagneticZ:real()                                 -- Imag Component                 imag_Ms_x = currents[c][pointCount].MagneticX:imag()                 imag_Ms_y = currents[c][pointCount].MagneticY:imag()                 imag_Ms_z = currents[c][pointCount].MagneticZ:imag()                 end                 -- Move to next vertice                 pointCount = pointCount+1             end                              -- Format average values as string              str = string.format("%15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G\n",                                                 currents.Axes[1][c],                                                 elementNum,                                                 centroidx/3,                                                 centroidy/3,                                                 centroidz/3,                                                 real_Js_x/3,                                                 real_Js_y/3,                                                 real_Js_z/3,                                                 imag_Js_x/3,                                                 imag_Js_y/3,                                                 imag_Js_z/3,                                                 real_Ms_x/3,                                                 real_Ms_y/3,                                                 real_Ms_z/3,                                                 imag_Ms_x/3,                                                 imag_Ms_y/3,                                                 imag_Ms_z/3)                                                  -- Print to file             file:write(str)                          -- Move to next triangle in solution             elementNum = elementNum + 1         end     end        end --Close the text file file:close() --return currents

 

Thanks!

Tagged:

Answers

  • Torben Voigt
    Torben Voigt
    Altair Employee
    edited August 2

    Hi Ricardo,

    Thanks for reporting! However, Feko 2020 is of course no longer up-to-date. Perhaps you have the opportunity to test it with Feko 2024?

    Best regards,
    Torben

  • Ricardo Sendrea
    Ricardo Sendrea Altair Community Member
    edited August 7
    Torben Voigt_20420 said:

    Hi Ricardo,

    Thanks for reporting! However, Feko 2020 is of course no longer up-to-date. Perhaps you have the opportunity to test it with Feko 2024?

    Best regards,
    Torben

    Hi Torben,

    Thanks for your suggestion!

    In the updated FEKO version, the .out currents match the POSTFEKO solution perfectly! I found a small bug in my export LUA script. Now my visualizations match FEKO exactly. For other's reference, I will include the fixed code here:

    app = pf.GetApplication() app:NewProject() app:OpenFile([[DIRECTORY\DESIGN.fek]]) -- Set up the file details absPath  = [[DIRECTORY\]]  -- HIGHLY recommended to hard code the absolute path! meshfile = "output.txt" --name of textile to use -- Find all possible data files names = pf.SurfaceCurrentsAndCharges.GetNames() --Get characteristic mode model currents currents = pf.SurfaceCurrentsAndCharges.GetDataSet(names[1]) --Load a given model sConf = app.Models["MODEL_NAME"].Configurations[1] --Import the mesh mesh = sConf.Mesh count = mesh.TriangleFaces.Count --Number of Triangles --Get mesh values points = mesh.Points --I/O for mesh write upper-- Open a file to write to. file = io.open(absPath..meshfile, "w") --Create Column Headers Vertice Index - Position(X,Y,Z) colHeader = string.format("%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s\n", "Frequency", "Element","X", "Y", "Z", "real_Jx", "realJy", "real_Jz", "imag_Jx", "imag_Jy", "imag_Jz", "real_Mx", "real_My", "real_Mz", "imag_Mx", "imag_My", "imag_Mz") --Write Column header to file file:write(colHeader) -- Get vector of all solved frequencies frequencies = currents.Axes[1] for c = 1, #frequencies do         str = ""         elementNum = 1         pointCount = 1         for i = 1, count do             --Load mesh triangle             triangle = mesh.TriangleFaces[i].Triangles             count2 = triangle.Count --Number of triangle vertices               for j = 1, count2 do             info = triangle[j].VertexIndices                                          -- Pre-Allocate Center of Triangle Value for ouput              -- Electric Current             real_Js_x = 0             real_Js_y = 0             real_Js_z = 0             imag_Js_x = 0             imag_Js_y = 0             imag_Js_z = 0                          -- Magnetic Current             real_Ms_x = 0             real_Ms_y = 0             real_Ms_z = 0             imag_Ms_x = 0             imag_Ms_y = 0             imag_Ms_z = 0                          -- Element Centroid             centroidx = 0             centroidy = 0             centroidz = 0                          for ii = 1, #info do                                                    -- Get vertice location                 position = points[triangle[j].VertexIndices[ii]]                 -- Add to centroid                 centroidx = centroidx + position[1]                 centroidy = centroidy + position[2]                 centroidz = centroidz + position[3]                                  -- Add Electric Currents                 real_Js_x = currents[c][pointCount].ElectricX:real() + real_Js_x                 real_Js_y = currents[c][pointCount].ElectricY:real() + real_Js_y                 real_Js_z = currents[c][pointCount].ElectricZ:real() + real_Js_z                                 -- Imag component                 imag_Js_x = currents[c][pointCount].ElectricX:imag() + imag_Js_x                 imag_Js_y = currents[c][pointCount].ElectricY:imag() + imag_Js_y                 imag_Js_z = currents[c][pointCount].ElectricZ:imag() + imag_Js_z                                      if currents[c][pointCount].MagneticX ~= nil then                 -- Add Magnetic Currents                 real_Ms_x = currents[c][pointCount].MagneticX:real() + real_Ms_x                 real_Ms_y = currents[c][pointCount].MagneticY:real() + real_Ms_y                 real_Ms_z = currents[c][pointCount].MagneticZ:real() + real_Ms_z                                 -- Imag Component                 imag_Ms_x = currents[c][pointCount].MagneticX:imag() + imag_Ms_x                 imag_Ms_y = currents[c][pointCount].MagneticY:imag() + imag_Ms_y                 imag_Ms_z = currents[c][pointCount].MagneticZ:imag() + imag_Ms_z                 end                 -- Move to next vertice                 pointCount = pointCount+1             end                              -- Format average values as string              str = string.format("%15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G\n",                                                 currents.Axes[1][c],                                                 elementNum,                                                 centroidx/3,                                                 centroidy/3,                                                 centroidz/3,                                                 real_Js_x/3,                                                 real_Js_y/3,                                                 real_Js_z/3,                                                 imag_Js_x/3,                                                 imag_Js_y/3,                                                 imag_Js_z/3,                                                 real_Ms_x/3,                                                 real_Ms_y/3,                                                 real_Ms_z/3,                                                 imag_Ms_x/3,                                                 imag_Ms_y/3,                                                 imag_Ms_z/3)                                                  -- Print to file             file:write(str)                          -- Move to next triangle in solution             elementNum = elementNum + 1         end     end        end --Close the text file file:close() --return currents
  • Torben Voigt
    Torben Voigt
    Altair Employee
    edited August 8
    Ricardo Sendrea said:

    Hi Torben,

    Thanks for your suggestion!

    In the updated FEKO version, the .out currents match the POSTFEKO solution perfectly! I found a small bug in my export LUA script. Now my visualizations match FEKO exactly. For other's reference, I will include the fixed code here:

    app = pf.GetApplication() app:NewProject() app:OpenFile([[DIRECTORY\DESIGN.fek]]) -- Set up the file details absPath  = [[DIRECTORY\]]  -- HIGHLY recommended to hard code the absolute path! meshfile = "output.txt" --name of textile to use -- Find all possible data files names = pf.SurfaceCurrentsAndCharges.GetNames() --Get characteristic mode model currents currents = pf.SurfaceCurrentsAndCharges.GetDataSet(names[1]) --Load a given model sConf = app.Models["MODEL_NAME"].Configurations[1] --Import the mesh mesh = sConf.Mesh count = mesh.TriangleFaces.Count --Number of Triangles --Get mesh values points = mesh.Points --I/O for mesh write upper-- Open a file to write to. file = io.open(absPath..meshfile, "w") --Create Column Headers Vertice Index - Position(X,Y,Z) colHeader = string.format("%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s\n", "Frequency", "Element","X", "Y", "Z", "real_Jx", "realJy", "real_Jz", "imag_Jx", "imag_Jy", "imag_Jz", "real_Mx", "real_My", "real_Mz", "imag_Mx", "imag_My", "imag_Mz") --Write Column header to file file:write(colHeader) -- Get vector of all solved frequencies frequencies = currents.Axes[1] for c = 1, #frequencies do         str = ""         elementNum = 1         pointCount = 1         for i = 1, count do             --Load mesh triangle             triangle = mesh.TriangleFaces[i].Triangles             count2 = triangle.Count --Number of triangle vertices               for j = 1, count2 do             info = triangle[j].VertexIndices                                          -- Pre-Allocate Center of Triangle Value for ouput              -- Electric Current             real_Js_x = 0             real_Js_y = 0             real_Js_z = 0             imag_Js_x = 0             imag_Js_y = 0             imag_Js_z = 0                          -- Magnetic Current             real_Ms_x = 0             real_Ms_y = 0             real_Ms_z = 0             imag_Ms_x = 0             imag_Ms_y = 0             imag_Ms_z = 0                          -- Element Centroid             centroidx = 0             centroidy = 0             centroidz = 0                          for ii = 1, #info do                                                    -- Get vertice location                 position = points[triangle[j].VertexIndices[ii]]                 -- Add to centroid                 centroidx = centroidx + position[1]                 centroidy = centroidy + position[2]                 centroidz = centroidz + position[3]                                  -- Add Electric Currents                 real_Js_x = currents[c][pointCount].ElectricX:real() + real_Js_x                 real_Js_y = currents[c][pointCount].ElectricY:real() + real_Js_y                 real_Js_z = currents[c][pointCount].ElectricZ:real() + real_Js_z                                 -- Imag component                 imag_Js_x = currents[c][pointCount].ElectricX:imag() + imag_Js_x                 imag_Js_y = currents[c][pointCount].ElectricY:imag() + imag_Js_y                 imag_Js_z = currents[c][pointCount].ElectricZ:imag() + imag_Js_z                                      if currents[c][pointCount].MagneticX ~= nil then                 -- Add Magnetic Currents                 real_Ms_x = currents[c][pointCount].MagneticX:real() + real_Ms_x                 real_Ms_y = currents[c][pointCount].MagneticY:real() + real_Ms_y                 real_Ms_z = currents[c][pointCount].MagneticZ:real() + real_Ms_z                                 -- Imag Component                 imag_Ms_x = currents[c][pointCount].MagneticX:imag() + imag_Ms_x                 imag_Ms_y = currents[c][pointCount].MagneticY:imag() + imag_Ms_y                 imag_Ms_z = currents[c][pointCount].MagneticZ:imag() + imag_Ms_z                 end                 -- Move to next vertice                 pointCount = pointCount+1             end                              -- Format average values as string              str = string.format("%15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G %15G\n",                                                 currents.Axes[1][c],                                                 elementNum,                                                 centroidx/3,                                                 centroidy/3,                                                 centroidz/3,                                                 real_Js_x/3,                                                 real_Js_y/3,                                                 real_Js_z/3,                                                 imag_Js_x/3,                                                 imag_Js_y/3,                                                 imag_Js_z/3,                                                 real_Ms_x/3,                                                 real_Ms_y/3,                                                 real_Ms_z/3,                                                 imag_Ms_x/3,                                                 imag_Ms_y/3,                                                 imag_Ms_z/3)                                                  -- Print to file             file:write(str)                          -- Move to next triangle in solution             elementNum = elementNum + 1         end     end        end --Close the text file file:close() --return currents

    Hi Ricardo,

    Thank you, I'm glad the problem has been fixed!

    Best regards,
    Torben

Categories