Problems/Confusion on Coding Bilinear and Linear Form

General
1

Screenshot 2024-05-23 022613
Screenshot 2024-05-23 022653
Hello,

I want to write the bilinear and linear forms to solve my problem.Desired bilinear and linear forms are attached to this message as an image. I have two main questions about implementing these equations:

  1. How can I defuned the unit normal vector and calculate the inner product between the unit normal vector and my Nedelec basis functions (for the third integral in the summation term of bilinear form)

  2. How can I define integrals over the facets. Ok, in the tutorials, I see generally they used ‘ds’ to calculate integrals because of the boundaries: ds = ufl.Measure(“ds”, domain, subdomain_data=facet_tags)

For example, for my case, my problem is 2D and integrals over the boundaries (input and output ports) are 1D. How can I define these integra, how should I write the code. For example, for bilenear form’s third integral, I will calculate integral of the defined inner product over the input and output ports. For 2D case, boundaries ar einput and putput edges of the rectangular waveguide. Also, for this example, two port’s effects will be inserted to calculations because integral is defined over the summation of integrals obtained from input and output ports, How can I define input and output ports and calculate the integrals correctly? Also, how can I defined these ports (these two ports (input and output) are different ports, how can I handle it) and calculate the inner products over selected ports? I am really confused about implementing these forms but I know it is possible in Fenicsx, I can implement these equations easily. Could you help me to wirte the code of whole two images, obtained equations as a A and b and the I will calculate the results, unknown E-fields.

Also, using these obtained bilinear and linear forms, I can calculate the my E fields, my aim, similarly with the Poisson equation example, right? from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, bcs=[bc], petsc_options={“ksp_type”: “preonly”, “pc_type”: “lu”})
uh = problem.solve()

Thank you very much

3

Given a mesh, you can get the normal to the facets with:

n = ufl.FacetNormal(mesh)

As for integrals over facets, you indeed need to use a ds measure. If your mesh is 2D then you can integrate over your 1D boundaries using a ds. You’d have to tag your input and output ports (for instance 1 for the input port and 2 for the output port), then you can integrate over the input and output with elements ds(1) and ds(2).

5

Thank you very much for your answer,

I implemented this problem (equations given in the above), I didn’t encounter any errors, but the results are weird, not in line with theoretical expectations. Do you see any obvious mistakes in the implementation?
I will add import and mesh generation steps additionally, and in the below of these import and my mesh generatio codes, I will give the main body of the code. I may have imported some packages more than once in the importing step, but this code works well with no error. I used Dolfinx 0.7.2 and Dolfinx Complex Python 3 Jupyter Notebook.

from mpi4py import MPI
from petsc4py import PETSc
import gmsh
import numpy as np
import ufl
from dolfinx.mesh import locate_entities_boundary, meshtags
from dolfinx.io import gmshio
from dolfinx import fem, mesh, plot
import pyvista
import scipy.sparse
import dolfinx
from dolfinx import mesh, fem, plot, io
import dolfinx.fem.petsc
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant,  Function, functionspace, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological,assemble_matrix)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.plot import vtk_mesh
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, 
                 div, dot, dx, grad, inner, lhs, rhs)
import scipy
from dolfinx import mesh as dfx_mesh
# Define geometric and topological dimensions
gdim = 2
shape = "triangle"
degree = 1

# Define the cell and domain
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))

# Define all the vertices for the mesh
vertices = np.array([[0.1, 0.1],  # Vertex 0
                     [0.1, 1.1],  # Vertex 1
                     [0.6, 0.1],  # Vertex 2
                     [0.6, 1.1],  # Vertex 3
                     [1.1, 0.1],  #  Vertex 4
                     [1.1, 1.1]])  # Vertex 5
vertices=vertices*(10**(-2))                  
# Define four cells (triangles) using the vertices defined above
cells = np.array([[0, 1, 2],   #  1 triangle
                  [1, 2, 3],   #  2 triangle
                  [2, 4, 3],   #  3 triangle
                  [4, 5, 3]],  #  4 triangle
                 dtype=np.int32)

# Create the mesh
meshc = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, vertices, domain)
tdim = meshc.topology.dim
fdim = tdim - 1
meshc.topology.create_connectivity(fdim, tdim)

# Check if PyVista backend is available
print(pyvista.global_theme.jupyter_backend)
pyvista.start_xvfb()  # Start virtual framebuffer if necessary

# Create and plot the mesh using PyVista
topology, cell_types, geometry = plot.vtk_mesh(meshc, gdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()

# Show the plot in an interactive window or save to a file
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("adjacent_triangles_mesh.png")
import gmsh
from dolfinx import mesh as dfx_mesh
from dolfinx.io import gmshio
import pyvista
from dolfinx import plot

def refine_mesh(cell_to_node, coords, num_refinements):
    """Refine the mesh recursively."""
    if num_refinements == 0:
        # Base case, return the original cell to node connectivity and coordinates
        return cell_to_node, coords
    else:
        # Create new triangles using midpoints
        new_cell_to_node = []
        new_coords = coords.copy()
        # Dictionary to keep track of the midpoints that have already been added
        edge_to_midpoint = {}

        for cell in cell_to_node:
            # Find the midpoints for the current cell
            cell_coords = coords[cell]
            midpoints = []
            midpoint_indices = []

            for i in range(len(cell)):
                edge = tuple(sorted((cell[i], cell[(i + 1) % len(cell)])))

                if edge not in edge_to_midpoint:
                    # Calculate midpoint for a new edge
                    midpoint = (cell_coords[i] + cell_coords[(i + 1) % len(cell)]) / 2.0
                    new_coords = np.vstack((new_coords, midpoint))
                    midpoint_index = len(new_coords) - 1
                    edge_to_midpoint[edge] = midpoint_index
                else:
                    # Get the existing midpoint index for the edge
                    midpoint_index = edge_to_midpoint[edge]

                midpoints.append(new_coords[midpoint_index])
                midpoint_indices.append(midpoint_index)

            # Define new cells using the original vertices and midpoints
            new_cells = [
                [cell[0], midpoint_indices[0], midpoint_indices[2]],
                [midpoint_indices[0], cell[1], midpoint_indices[1]],
                [midpoint_indices[2], midpoint_indices[1], cell[2]],
                midpoint_indices,  # the central triangle
            ]
            new_cell_to_node.extend(new_cells)

        # Recursively refine the new mesh
        return refine_mesh(new_cell_to_node, new_coords, num_refinements - 1)

# Continue with the rest of your code to use the refined mesh...


# Get the cell to node connectivity and coordinates of the original mesh
cell_to_node = meshc.topology.connectivity(2, 0).array.reshape(-1, 3)
coords = meshc.geometry.x

# Choose the number of refinements
num_refinements = 1  # Change this to the desired number of refinement levels

# Refine the mesh
refined_cell_to_node, refined_coords = refine_mesh(cell_to_node, coords, num_refinements)

cells = np.asarray(refined_cell_to_node, dtype=np.int32)
xs = np.asarray(refined_coords, dtype=np.float64).reshape(-1, 3)  
pos = refined_coords[:, :2]
# Define geometric and topological dimensions
gdim = 2
shape = "triangle"
degree = 1

# Define the cell and domain
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
# Create the mesh
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, refined_cell_to_node, pos, domain)
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)

# Check if PyVista backend is available
print(pyvista.global_theme.jupyter_backend)
pyvista.start_xvfb()  # Start virtual framebuffer if necessary

# Create and plot the mesh using PyVista
topology, cell_types, geometry = plot.vtk_mesh(mesh, gdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()

# Show the plot in an interactive window or save to a file
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("adjacent_triangles_mesh.png")

######MAIN BODY OF THE CODE/IMPLEMENTATION 

#Spatial Coordinates, Function Spaces and Facet Normal

x = SpatialCoordinate(mesh)
V = FunctionSpace(mesh, ("Nedelec 1st kind H(curl)", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
n = ufl.FacetNormal(mesh)


# %%%WR 975 WAVEGUIDE Parameters%%%
length_a_waveguide= (247.65)*(10**(-3))
length_b_waveguide= (123.825)*(10**(-3))
frequency=1*(10**9)
speed_of_light=(3)*(10**8)
wavelength=(speed_of_light)/frequency
k0=(2*(np.pi))/wavelength
# angular_freq=2*np.pi*frequency
# first_term=((angular_freq)**2)/((speed_of_light)**2)
first_term=k0**2
second_term= ((np.pi)/length_a_waveguide)**2
beta_10=first_term-second_term
beta_10

# ##Start to weak form
L_lhs= (ufl.inner(ufl.curl(u),ufl.curl(v))*ufl.dx) - ((beta_10)*(beta_10)*(ufl.inner(u,v)*ufl.dx)) 

coordinates = mesh.geometry.x[:, :2] 
# Find the minimum and maximum values for x and y
min_x, min_y = np.min(coordinates, axis=0)
max_x, max_y = np.max(coordinates, axis=0)

# Create boundary conditions using these values
boundaries = [
    (1, lambda x: np.isclose(x[0], min_x)),
    (2, lambda x: np.isclose(x[0], max_x)),
    (3, lambda x: np.isclose(x[1], min_y)),
    (4, lambda x: np.isclose(x[1], max_y))
]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

#Dirichlet BCs

facets3 = facet_tag.find(3)
u_D3 = fem.Function(V)
u_D3.x.array[:] = 0
dofs_D3 = fem.locate_dofs_topological(V, fdim, facets3)
bcD3 = fem.dirichletbc(u_D3, dofs_D3)

facets4 = facet_tag.find(4)
u_D4 = fem.Function(V)
u_D4.x.array[:] = 0
dofs_D4 = fem.locate_dofs_topological(V, fdim, facets4)
bcD4 = fem.dirichletbc(u_D4, dofs_D4)

bc = [bcD3, bcD4]

# Definition of cross product in 2D
n = ufl.FacetNormal(mesh)
def cross_2D(n, v):
    #  n = (0, 0, 1) for a mesh in the xy-plane
    return ufl.as_vector((-v[1], v[0]))


# Input port
bc1 = 1j * (beta_10) * ufl.inner(cross_2D(n, u), cross_2D(n,v)) * ds(1)
L_lhs += bc1
# Output port
bc2 = 1j * (beta_10) * ufl.inner(cross_2D(n, u), cross_2D(n,v)) * ds(2)
L_lhs += bc2


#E_inc definition

class Einc:
    def __init__(self, E0, a):
        self.E0 = E0
        self.a = a
    def eval(self, x):
        ey = self.E0 * np.sin(np.pi * x[0] / self.a)
        ex = np.full_like(ey, 0.0+0j, dtype=np.complex128)
        return (ex,ey)

E0 = 0.5
incident_field = Einc(E0, length_a_waveguide)

#Build RHS
uh = fem.Function(V, dtype=np.complex128)
uh.interpolate(incident_field.eval)
f_inc = 2j * (beta_10) * ufl.inner(cross_2D(n, u), cross_2D(n, uh)) * ds(1)
uh.name = "Einc"

##Solving the equation
problem = fem.petsc.LinearProblem(L_lhs, f_inc, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
solution=problem.solve()
print(solution)
6

Please revise your code and make aninimal reproducible example. There are so many imports here that are duplicate.

Please also remove all code that is commented out, and maybe add a visualization of the solution you get, and some comments about what you expect.

7

image
image935×698 63.4 KB

E_h_projected = fem.Function(W)
E_h_magnitude_expr = ufl.sqrt(ufl.inner(solution, solution))
E_h_projected.interpolate(fem.Expression(E_h_magnitude_expr, W.element.interpolation_points(), mesh.comm))

# Create the PyVista grid
topology, cell_types, geometry = plot.vtk_mesh(mesh, mesh.geometry.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.point_data["E_h"] = E_h_projected.vector.array

# Plot the solution
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, scalars="E_h", show_scalar_bar=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    plotter.screenshot("solution.png")

I edited my previous message, I removed duplicated imports. Then, I visualized my solution with this code. When I visualize the result (image is attached to this message), I would expect the result here to be symmetrical with respect to the x-axis. In other words, there is a more energy density at the end of the x-axis, and a less at the input port.

8

Please provide a minimum working example. Your previous code is very long, with many redundant imports/code lines and many lines of code irrelevant to the issue, like mesh refinement. You may also want to simply create your mesh using dolfinx.mesh.create_rectangle rather than manually. This will make it easy for other people to run your code on their end and help you.