Paraview visualization

I/O
1

Hello all,

I am trying to visualize some eigenmodes I calculated with fenics.
I work with gmsh to generate my meshes:

import gmsh
import sys

lc=0.08
gmsh.initialize()
gmsh.model.add("Nanowire")
L,B,H,r = 1,1,1,0.2
gmsh.model.occ.addBox(0, 0, 0, L, B, H,1)
gmsh.model.occ.addCylinder(0.5, 0.5,0.1,0, 0, 0.8 , r, 2)
gmsh.model.occ.synchronize()

#Refine cylinder mesh
gmsh.model.mesh.field.add("Cylinder", 3)
gmsh.model.mesh.field.setNumber(3, "VIn", lc / 3)
gmsh.model.mesh.field.setNumber(3, "VOut", lc)
gmsh.model.mesh.field.setNumber(3, "Radius", 0.2)
gmsh.model.mesh.field.setNumber(3, "XCenter", 0.5)
gmsh.model.mesh.field.setNumber(3, "YCenter", 0.5)
gmsh.model.mesh.field.setNumber(3, "ZCenter", 0.1)
gmsh.model.mesh.field.setNumber(3, "XAxis", 0)
gmsh.model.mesh.field.setNumber(3, "YAxis", 0)
gmsh.model.mesh.field.setNumber(3, "ZAxis", 0.8)

gmsh.model.mesh.field.setAsBackgroundMesh(3)

gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)



gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(3)

gmsh.write("Nanowire.msh")

if '-nopopup' not in sys.argv:
    gmsh.fltk.run()
    
gmsh.finalize()

I convert it to xdmf via

import meshio


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:geometrical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh


msh = meshio.read("Nanowire.msh")
print(msh)

#triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
#meshio.write("Nanowire2D_facets.xdmf", triangle_mesh)

tetra_mesh = create_mesh(msh, "tetra", prune_z=False)
meshio.write("Nanowire_tetra.xdmf", tetra_mesh)

and the use the mesh as

    import numpy as np
    import matplotlib.pyplot as plt

    N=40
    mesh = UnitCubeMesh(N,N,N)
    with XDMFFile("Nanowire_tetra.xdmf") as file:
        file.read(mesh)
    V = VectorFunctionSpace(mesh, family='CG', degree=1, dim=3)
    bc = DirichletBC(V,[0, 0, 0],boundary)

    materials = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
    eps = Permittivity(materials, 4, 1, degree=0)
    subdomain_0 = Omega_0()
    subdomain_1 = Omega_1()
    subdomain_0.mark(materials, 0)
    subdomain_1.mark(materials, 1)

    space = helmholtz_wellposed(mesh, V, bc, 5E3, eps)
    print(len(space))
    for i,eigfun in enumerate(space):
        vtkfile = File(str(i)+'.pvd')
        vtkfile << eigfun

The problem statement is as follows:

    u, v = TrialFunction(V), TestFunction(V)
    a = inner(curl(u), curl(v))*dx(domain=mesh)
    L = inner(Constant((0, 0, 0)), v)*dx(domain=mesh)
    m = eps*inner(u,v)*dx(domain=mesh)
    A, _ = assemble_system(a, L, bc)
    B = assemble(m)

I am simulating a cylinder (lets call it nanowire) and I do get modes that conform to the symmetry, but there is some weird speckle-like pattern in the ParaView visualization, which can be seen here (a cut):

s1
s11074×851 41.8 KB

Does anybody know what this pattern is and how I can get rid of it? It changes shape when I turn the camera, so it is likely something to do with Paraview.

I am running linux:

$uname -v
#1 SMP Debian 4.19.235-1 (2022-03-17)
$ sudo lshw -C video
  *-display                 
       description: VGA compatible controller
       product: Intel Corporation
       vendor: Intel Corporation
       physical id: 2
       bus info: pci@0000:00:02.0
       version: 00
       width: 64 bits
       clock: 33MHz
       capabilities: pciexpress msi pm vga_controller bus_master cap_list rom
       configuration: driver=i915 latency=0
       resources: irq:128 memory:a1000000-a1ffffff memory:b0000000-bfffffff ioport:6000(size=64) memory:c0000-dffff
  *-display UNCLAIMED
       description: Display controller
       product: Topaz XT [Radeon R7 M260/M265 / M340/M360 / M440/M445]
       vendor: Advanced Micro Devices, Inc. [AMD/ATI]
       physical id: 0
       bus info: pci@0000:01:00.0
       version: c3
       width: 64 bits
       clock: 33MHz
       capabilities: pm pciexpress msi cap_list
       configuration: latency=0
       resources: memory:90000000-9fffffff memory:a0000000-a01fffff ioport:5000(size=256) memory:a2300000-a233ffff memory:a2340000-a235ffff
2

Looks like when paraview tries to plot overlapping data. Are you sure you only have one data block in your pipeline?

3

I think so, find attached a screenshot of the whole GUI:

Bildschirmfoto von 2022-05-28 08-53-26
Bildschirmfoto von 2022-05-28 08-53-261920×1080 198 KB

4

I think it might not be an issue with visualization, but with my function space. Trying to change the definition to Nedelec-Polynomials, I get the following error:

*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  Illegal value rank (1), expecting (2).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
5

I would suggest you make a minimal example that can reproduce the error, as it is hard to give guidance with only snippets of code.

6

Hi dokken,

Gladly. This is my code. I guess most of it might look familiar to you, as I copied it from various parts of the documentation :wink:

from dolfin import *
from dolfin.cpp.io import XDMFFile

def orthogonalize(A):
    """L^2-orthogonalize a list of Functions living on the same
    function space. Modify the functions in-place.
    Use classical Gramm-Schmidt algorithm for brevity.
    For numerical stability modified Gramm-Schmidt would be better.
    """

    # Set of single function is orthogonal
    if len(A) <= 1:
        return

    # Orthogonalize overything but the last function
    orthogonalize(A[:-1])

    # Orthogonalize the last function to the previous ones
    f = A[-1]
    for v in A[:-1]:
        r = assemble(inner(f, v)*dx) / assemble(inner(v, v)*dx)
        assert f.function_space() == v.function_space()
        f.vector().axpy(-r, v.vector())



def helmholtz_wellposed(mesh, V, bc, approxEig, eps):
    u, v = TrialFunction(V), TestFunction(V)
    a = inner(curl(u), curl(v))*dx(domain=mesh)
    L = inner(Constant((0, 0, 0)), v)*dx(domain=mesh)
    m = eps*inner(u,v)*dx(domain=mesh)
    A, _ = assemble_system(a, L, bc)
    B = assemble(m)

    eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B))
    eigensolver.parameters['problem_type'] = 'gen_hermitian'
    eigensolver.parameters['spectrum'] = 'target real'
    eigensolver.parameters['spectral_shift'] = approxEig
    eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
    eigensolver.parameters['tolerance'] = 1e-6
    #eigensolver.parameters['verbose'] = True
    eigensolver.solve(10) # Find 5 eigenpairs close to target
    space = []
    for j in range(eigensolver.get_number_converged()):
        r, c, rx, cx = eigensolver.get_eigenpair(j)
        if near(r, approxEig, approxEig/2):
            print(j,r)
            eig = Function(V)
            eig.vector()[:] = rx
        space.append(eig)

    orthogonalize(space)
    return space


def boundary(x):
    tol = 1E-14
    return near(x[0], 0, tol) or near(x[1], 0, tol) or near(x[2], 0, tol) \
        or near(x[0], 1, tol) or near(x[1], 1, tol) or near(x[2], 1, tol) \

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        r = np.sqrt((x[0]-0.5)**2+(x[1]-0.5)**2)
        return r <= 0.2 + tol

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        r = np.sqrt((x[0]-0.5)**2+(x[1]-0.5)**2)
        return r >= 0.2 - tol

class Permittivity(UserExpression):
    def __init__(self, materials, f_0, f_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.f_0 = f_0
        self.f_1 = f_1

    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.f_0
        else:
            values[0] = self.f_1



if __name__ == '__main__':
    import numpy as np
    import matplotlib.pyplot as plt

    mesh = Mesh()
    with XDMFFile("Nanowire_tetra.xdmf") as file:
        file.read(mesh)
    #V = VectorFunctionSpace(mesh, family='CG', degree=1, dim=3)
    V = VectorFunctionSpace(mesh, family='RT', degree=2)
    bc = DirichletBC(V,[0, 0, 0], boundary)

    materials = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
    eps = Permittivity(materials, 4, 1, degree=0)
    subdomain_0 = Omega_0()
    subdomain_1 = Omega_1()
    subdomain_0.mark(materials, 0)
    subdomain_1.mark(materials, 1)

    space = helmholtz_wellposed(mesh, V, bc, 5E3, eps)
    print(len(space))
    for i,eigfun in enumerate(space):
        vtkfile = File(str(i)+'.pvd')
        vtkfile << eigfun
7

I managed to use the N1curl elements by defining a function space via

V = FunctionSpace(mesh, 'N1curl', 3)

The boundary conditions stated above are accepted.
The solution still had the pattern I showed above; I am starting to suspect the mesh to be the issue (similarly as in other cases in this forum where poeple solved MWEQ, although slightly differently to what I am doing) and reverted to a UnitCubeMesh. I see what comes out.

8

Note that using VTKFile will interpolate the solution into CG-1, which could explain your visualization issue.
I would recommend:

  1. Interpolate solution into a suitable vector DG space
  2. Use XDMFFile.write_checkpoint to save the solution.
9

Thank you for the suggestion, but replacing the vtk output with these lines:

    V_plot = VectorFunctionSpace(mesh, 'DG',degree=1,dim=3)

    for i,eigfun in enumerate(space):
        eig_int = interpolate(eigfun,V_plot)
        with XDMFFile(str(i)+".xdmf") as xdmf:
            xdmf.write(eig_int)

did not help. I also tried out RT elements.

10

I said

Not xdmf.write as this function also interpolates into CG 1

11

Going to write_checkpoint does not improve things.
This speckle-like pattern persists for the non-uniform mesh. It disappears for the uniform mesh, but even there in the solution has some problems:

  1. The discretization is clearly visible, even when using RT3 elements and interpolating down to DG1 (The number should indicate the degree, I hope this notation is correct):
    V = FunctionSpace(mesh, 'RT', 3)
    bc = DirichletBC(V,[0, 0, 0], boundary)

    materials = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
    eps = Permittivity(materials, 4, 1, degree=0)
    subdomain_0 = Omega_0()
    subdomain_1 = Omega_1()
    subdomain_0.mark(materials, 0)
    subdomain_1.mark(materials, 1)

    space = helmholtz_wellposed(mesh, V, bc, 1E14, eps)

    V_plot = VectorFunctionSpace(mesh, 'DG',degree=1,dim=3)

    for i,eigfun in enumerate(space):
        eig_int = interpolate(eigfun,V_plot)
        with XDMFFile(str(i)+".xdmf") as xdmf:
            xdmf.write_checkpoint(eig_int, "sol")
  1. The solution seems to be constant over the cylinder volume- a waveguide mode should show a maximum in the center of the wire and then decrease towards the wire boundaries like a gaussian.

I do not know, at this point, what the issues stated above is caused by. It could be the problem statement. What is clear to me now is that the pattern I am talking about disappears for the uniform mesh. I could live with the speckle pattern if the problems above are somehow solved. I will open a new discussion topic detailing the steps I took for deriving at this problem statement.