Inaccurate interpolation of function

General
1

Hi,

I am modeling a 3D magnetostatic problem in which I want to describe the current density in a circular coil using Nedelec elements. The current density is circular and I define it myself. When I visualize the it by interpolating from the original Nedelec to DG space, I see something very irregular. What is more worrisome, even the Jz component is very strongly-nonzero at points when it is originally defined to be identically zero everywhere. Below is a screenshot of the paraview visualization of the magnitude distribution (min, max should be 0,1 ideally):

image
image964×912 35.3 KB

The function I use to define the current density gives the ideal distribution whose magnitude looks like this:

image

I am unsure where the problem lies, if at all. Is this is a bug/inaccuracy in the interpolation process in dolfinx, how things can be best fit into these mathematical functions or Paraview is doing a poor job at visualization?

Below is an MWE that generates the Js.xdmf file which can be visualized in Paraview. It should be run in dolfinx 0.6.0 because the latest release somehow does not work with the last part of saving as xdmf file

import dolfinx, gmsh, mpi4py, ufl
import numpy as np
import matplotlib.pyplot as plt
from dolfinx.io import gmshio

cyl_rad = 0.1, 0.2
cyl_height = 0.25
fov_rad = 0.5

#%% create a cylindrical shell inside a box
def create_geom(disp=False):
    gmsh.initialize()
    gmsh.clear()
    gdim = 3
    occ = gmsh.model.occ
    # inner cylinder
    obj_cyl_disk = occ.addDisk(0, 0, 0, cyl_rad[0], cyl_rad[0])
    _, obj_cyl1, _ = occ.extrude([(gdim-1, obj_cyl_disk)], 0, 0, cyl_height)
    # outer cylinder
    obj_cyl_disk = occ.addDisk(0, 0, 0, cyl_rad[1], cyl_rad[1])
    _, obj_cyl2, _ = occ.extrude([(gdim-1, obj_cyl_disk)], 0, 0, cyl_height)
    # enclosure
    obj_fov = occ.addBox(-fov_rad, -fov_rad, 0, 2*fov_rad, 2*fov_rad, fov_rad)
    final = occ.fragment([(gdim, obj_fov)], [obj_cyl1, obj_cyl2])

    occ.synchronize()
    # tag domains and boundaries
    all_doms = gmsh.model.getEntities(gdim)
    for idx, dom in enumerate(all_doms):
        gmsh.model.addPhysicalGroup(dom[0], (dom[1],), idx+1)
    all_edges = gmsh.model.getEntities(gdim-1)
    for idx, edge in enumerate(all_edges):
        gmsh.model.addPhysicalGroup(edge[0], (edge[1],), edge[1])

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.refine()

    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, 0)
    if disp:
        gmsh.fltk.run()
    gmsh.finalize()
    return mesh, ct, ft

mesh, ct, ft = create_geom(disp=False) # make True for visual inspection

# We have a circular coil, define a current density vector inside it
def current_density(r):
    '''define circular current density in the coil which is zero everywhere outside'''
    R_in, R_out = cyl_rad # read the global variable cyl_rad
    height = cyl_height # height, read from the global variable
    x, y, z = r
    rho = np.sqrt(x**2 + y**2)
    print("Interpolated coordinate radius min value = {}, max value = {}".format(rho.min(), rho.max()))
    print("Interpolated coordinates' min, max z value = {}, {}".format(z.min(), z.max()))

    theta = np.arctan2(y, x) # angle theta, relevant only in the rounded regions

    Jx = np.full_like(x, 0.0)
    Jy = np.full_like(x, 0.0)    
    Jz = np.full_like(x, 0.0) # this component will remain zero

    points_in_coil = (rho >= R_in)*(rho <= R_out)*(z>=0)*(z<=height)
    Jx[points_in_coil] = -np.sin(theta[points_in_coil])
    Jy[points_in_coil] = np.cos(theta[points_in_coil])

    return Jx, Jy, Jz

# for electromagnetic problem, we will need current density on curl elements
curl_el = ufl.FiniteElement("N1curl", mesh.ufl_cell(), 2)
Omega = dolfinx.fem.FunctionSpace(mesh, curl_el)
Js = dolfinx.fem.Function(Omega)
Js.interpolate(current_density)

# For visualization, interpolate to DG space
Omega_dg = dolfinx.fem.VectorFunctionSpace(mesh, ("DG", 2))
Js_dg = dolfinx.fem.Function(Omega)
Js_dg.interpolate(Js)

# save as xdmf to visualize in paraview
with dolfinx.io.XDMFFile(mesh.comm, "Js.xdmf", "w", encoding=dolfinx.io.XDMFFile.Encoding.ASCII) as file:
    file.write_mesh(mesh)
    file.write_function(Js_dg)

Thanks a lot for comments and suggestions.

1 Like
2

You should not use XDMFFile to visualize Nedelec functions.
XDMFile only supports (continuous) Lagrange. In your case, with 0.6.0 it interpolates the function into a first order continuous space.

You should interpolate the solution into a suitable DG space, and use VTXWriter or VTKFile to store the solution.

3

Thank you for your feedback!

I interpolated to DG but I guess xdmf is still not a suitable way to save. With VTXWriter, what is the correct reader in Paraview? It shows a huge list and so far everything I tried led to errors.

4

See Write mixed element function to files - #11 by dokken

5

Thank you so much!

I also have something to add for Windows users:

  • Paraview gives ADIOS2VTXReader option only in the MPI version. For some reason, the non-MPI version that I was using does not have it.
  • Secondly, if the is an error due to missing msmpi.dll, install MS MPI before installing Paraview’s MPI version.
1 Like