How to locate edges of a rectangle for defining dofs?

vishalkr · November 5, 2024, 7:09am

I have defined rectangle like this

upper_plate = gmsh.model.occ.addRectangle(2, 3, 0, 1, 0.1)
lower_plate = gmsh.model.occ.addRectangle(2, 2, 0, 1, 0.1)
gmsh.model.occ.synchronize()

V = fem.functionspace(mesh_obj, ("CG", 2))

# Define markers for plates
def upper_plate_boundary(x):
    return np.logical_and(np.logical_and((2 <= x[0]), (x[0]<=3)), np.logical_and((3 <= x[1]), (x[1]<=3.1))).any()

def lower_plate_boundary(x):
    return np.logical_and(np.logical_and((2 <= x[0]), (x[0]<=3)), np.logical_and((2 <= x[1]), (x[1]<=2.1))).any()



fdim = mesh_obj.topology.dim - 1
# Locate DOFs for the plates
upper_facets = locate_entities_boundary(mesh_obj, fdim, upper_plate_boundary)
lower_facets = locate_entities_boundary(mesh_obj, fdim, lower_plate_boundary)

But I am not able to locate edges of the rectangle to locate_entities_boundary. I also tried using cell tags, and face tags. But no sucess.

dokken · November 5, 2024, 7:11am

Please make a minimal reproducible example, as you clearly have attempted to solve this.

Giving only two snippets of your code increases the amount of effort anyone wanting to help you would have to make.

vishalkr · November 5, 2024, 7:25am

Here the code details

import numpy as np
import ufl
import gmsh
import dolfinx
from petsc4py import PETSc
from dolfinx import fem, mesh, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from dolfinx.io import gmshio, XDMFFile
from ufl import div, grad, inner, dx, Measure
from dolfinx.mesh import locate_entities_boundary, meshtags

# Initialize parameters
L, H = 5, 5  # Domain size
plate_length = 1
plate_height = 0.1
V0 = 100.0  # Voltage on upper plate
model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 2

# Initialize gmsh
gmsh.initialize()

if mesh_comm.rank == model_rank:
    # Create background domain
    background = gmsh.model.occ.addRectangle(0, 0, 0, L, H)
    
    # Create plates
    upper_plate = gmsh.model.occ.addRectangle(2, 3, 0, plate_length, plate_height)
    lower_plate = gmsh.model.occ.addRectangle(2, 2, 0, plate_length, plate_height)
    gmsh.model.occ.synchronize()
    
    # Fragment all geometries together
    all_surfaces = [(2, upper_plate), (2, lower_plate)]
    whole_domain = gmsh.model.occ.fragment([(2, background)], all_surfaces)
    gmsh.model.occ.synchronize()
    
    # Create physical markers
    background_surfaces = []
    plate_surfaces = []
    
    for domain in whole_domain[0]:
        com = gmsh.model.occ.getCenterOfMass(domain[0], domain[1])
        
        # Identify upper plate
        if 2 <= com[0] <= 3 and 3 <=com[1] <= 3.1:
            gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=1)
            gmsh.model.setPhysicalName(domain[0], 1, "upper_plate")
            plate_surfaces.append(domain)
            
        # Identify lower plate
        elif 2 <= com[0] <= 3 and 2 <=com[1] <= 2.1:
            gmsh.model.addPhysicalGroup(domain[0], [domain[1]], tag=2)
            gmsh.model.setPhysicalName(domain[0], 2, "lower_plate")
            plate_surfaces.append(domain)
            
        # Background is everything else
        else:
            background_surfaces.append(domain[1])
    
    # Add marker for background
    gmsh.model.addPhysicalGroup(2, background_surfaces, tag=0)
    gmsh.model.setPhysicalName(2, 0, "background")
    
    # Create refined mesh around plates
    gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.05)
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")
    gmsh.write("capacitor.msh")

# Read mesh and set up connectivity
mesh_obj, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
mesh_obj.topology.create_connectivity(gdim - 1, gdim)  # Create necessary connectivity
gmsh.finalize()

# Save mesh for visualization
with XDMFFile(MPI.COMM_WORLD, "capacitor.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_obj)
    xdmf.write_meshtags(cell_tags, mesh_obj.geometry)

# Create function space
V = fem.functionspace(mesh_obj, ("CG", 2))

# Define markers for plates
def upper_plate_boundary(x):
    return np.logical_and(np.logical_and((2 <= x[0]), (x[0]<=3)), np.logical_and((3 <= x[1]), (x[1]<=3.1))).any()

def lower_plate_boundary(x):
    return np.logical_and(np.logical_and((2 <= x[0]), (x[0]<=3)), np.logical_and((2 <= x[1]), (x[1]<=2.1))).any()



fdim = mesh_obj.topology.dim - 1
# Locate DOFs for the plates
upper_facets = locate_entities_boundary(mesh_obj, fdim, upper_plate_boundary)
lower_facets = locate_entities_boundary(mesh_obj, fdim, lower_plate_boundary)

# Create meshtags for the plates
upper_facets_tag = meshtags(mesh_obj, gdim, upper_facets, np.full_like(upper_facets, 1))
lower_facets_tag = meshtags(mesh_obj, gdim, lower_facets, np.full_like(lower_facets, 2))

# Locate DOFs based on the markers
upper_dofs = fem.locate_dofs_topological(V, gdim-1, upper_facets)
lower_dofs = fem.locate_dofs_topological(V, gdim-1, lower_facets)

# Apply Dirichlet boundary conditions
bc_upper = fem.dirichletbc(PETSc.ScalarType(V0), upper_dofs, V)
bc_lower = fem.dirichletbc(PETSc.ScalarType(-V0), lower_dofs, V)

# Combine boundary conditions
bcs = [bc_upper, bc_lower]

For mesh generation, I have followed ideas given in electromagnetic example of the tutorial - Electromagnetics example — FEniCSx tutorial.

I hope it is fine.

dokken · November 5, 2024, 9:04am

vishalkr:

def upper_plate_boundary(x):
    return np.logical_and(np.logical_and((2 <= x[0]), (x[0]<=3)), np.logical_and((3 <= x[1]), (x[1]<=3.1))).any()

def lower_plate_boundary(x):
    return np.logical_and(np.logical_and((2 <= x[0]), (x[0]<=3)), np.logical_and((2 <= x[1]), (x[1]<=2.1))).any()

First of all, remove the .any() on either call.

Secondly, note that an interior boundary is not a “boundary” in the sense of locate_entities_boundary.

I guess what you want to achieve is tagging the interface of what is marked between 0 and 1 and 0 and 2?

This is something that you can easily do through gmsh. See for instance: Meshes from external sources — FEniCS Tutorial @ Sorbonne

vishalkr · November 5, 2024, 6:36pm

yes, I understand I am selecting the area rather than boundaries.
Tried the suggested tutorial Meshes from external sources — FEniCS Tutorial @ Sorbonne

import numpy as np
import ufl
import gmsh
import dolfinx
from petsc4py import PETSc
from dolfinx import fem, mesh, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from dolfinx.io import gmshio, XDMFFile
from ufl import div, grad, inner, dx, Measure
from dolfinx.mesh import locate_entities_boundary, meshtags

# Initialize parameters
L, H = 5, 5  # Domain size
plate_length = 1
plate_height = 0.1
V0 = 100.0  # Voltage on upper plate
model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 2

# Initialize gmsh
gmsh.initialize()

if mesh_comm.rank == model_rank:
    # Create background domain
    background = gmsh.model.occ.addRectangle(0, 0, 0, L, H) #1
    
    # Create plates
    upper_plate = gmsh.model.occ.addRectangle(2, 3, 0, plate_length, plate_height) #2
    lower_plate = gmsh.model.occ.addRectangle(2, 2, 0, plate_length, plate_height) #3
    gmsh.model.occ.synchronize()
    
    # Fragment all geometries together
    all_surfaces = [(2, upper_plate), (2, lower_plate)]
    whole_domain, map_to_input = gmsh.model.occ.fragment([(2, background)], all_surfaces)
    gmsh.model.occ.synchronize()
    
    # Create physical markers
    background_surfaces = []
    plate_surfaces = []
    
    upper_mark = [idx for (dim, idx) in map_to_input[1] if dim == 2]
    lower_mark = [idx for (dim, idx) in map_to_input[2] if dim == 2]
    total_mark = upper_mark + lower_mark
    background_mark = [idx for (dim, idx) in map_to_input[0] if dim == 2 and idx not in total_mark]
    
    gmsh.model.addPhysicalGroup(2, upper_mark, tag=3)
    gmsh.model.addPhysicalGroup(2, lower_mark, tag=7)
    
    upper_boundary = gmsh.model.getBoundary([(2, e) for e in upper_mark], recursive=False, oriented=False)
    lower_boundary = gmsh.model.getBoundary([(2, e) for e in lower_mark], recursive=False, oriented=False)
    outer_boundary = gmsh.model.getBoundary([(2, e) for e in background_mark], recursive=False, oriented=False)
    
    upper_interface = [idx for (dim, idx) in upper_boundary if dim == 1]
    lower_interface = [idx for (dim, idx) in lower_boundary if dim == 1]
    total_interface = upper_interface + lower_interface
    ext_boundary = [idx for (dim, idx) in outer_boundary if idx not in total_interface and dim == 1]
    
    gmsh.model.addPhysicalGroup(1, upper_interface, tag=12)
    gmsh.model.addPhysicalGroup(1, lower_interface, tag=15)
    gmsh.model.addPhysicalGroup(1, ext_boundary, tag=17)
    
    # Create refined mesh around plates
    gmsh.model.mesh.field.add("Distance", 1)
    edges = gmsh.model.getBoundary(all_surfaces, oriented=False)
    gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
    
    # Add threshold field for mesh refinement
    gmsh.model.mesh.field.add("Threshold", 2)
    gmsh.model.mesh.field.setNumber(2, "IField", 1)
    gmsh.model.mesh.field.setNumber(2, "LcMin", plate_height/3)
    gmsh.model.mesh.field.setNumber(2, "LcMax", plate_height*6)
    gmsh.model.mesh.field.setNumber(2, "DistMin", plate_height*4)
    gmsh.model.mesh.field.setNumber(2, "DistMax", plate_height*10)
    
    gmsh.model.mesh.field.add("Min", 5)
    gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
    gmsh.model.mesh.field.setAsBackgroundMesh(5)
    
    # Generate mesh
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(2)
    gmsh.model.mesh.optimize("Netgen")
    gmsh.write("capacitor.msh")

# Read mesh and set up connectivity
mesh_obj, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, 0, gdim=2)
# mesh_obj.topology.create_connectivity(gdim - 1, gdim)  # Create necessary connectivity
gmsh.finalize()


# Save mesh for visualization
with XDMFFile(MPI.COMM_WORLD, "capacitor.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_obj)
    xdmf.write_meshtags(cell_tags, mesh_obj.geometry)

But getting error

MPI_ERR_RANK: invalid rank
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI COMMUNICATOR 9 DUP FROM 4
with errorcode 6.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------

Process ended with exit code 6.

vishalkr · November 6, 2024, 7:10am

Solved it, the problem was because physical group for background was not defined.

gmsh.model.addPhysicalGroup(2, background_mark, tag=9)

vishalkr · November 6, 2024, 11:56am

import numpy as np
import ufl
import gmsh
import dolfinx
from petsc4py import PETSc
from dolfinx import fem, mesh, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from dolfinx.io import gmshio, XDMFFile
from ufl import div, grad, inner, dx, Measure
from dolfinx.mesh import locate_entities_boundary, meshtags

# Initialize parameters
L, H = 3, 3  # Domain size
plate_length = 1
plate_height = 0.1
V0 = 1000.0  # Voltage on upper plate
model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 2



# Initialize gmsh
gmsh.initialize()

# Create background domain
background = gmsh.model.occ.addRectangle(0, 0, 0, L, H) #1

# Create plates
upper_plate = gmsh.model.occ.addRectangle(1, 2, 0, plate_length, plate_height) #2
lower_plate = gmsh.model.occ.addRectangle(1, 1.1, 0, plate_length, plate_height) #3
gmsh.model.occ.synchronize()

# Fragment all geometries together
all_surfaces = [(2, upper_plate), (2, lower_plate)]
whole_domain, map_to_input = gmsh.model.occ.fragment([(2, background)], all_surfaces)
gmsh.model.occ.synchronize()

# Create physical markers
background_surfaces = []
plate_surfaces = []

upper_mark = [idx for (dim, idx) in map_to_input[1] if dim == 2]
lower_mark = [idx for (dim, idx) in map_to_input[2] if dim == 2]
total_mark = upper_mark + lower_mark
background_mark = [idx for (dim, idx) in map_to_input[0] if dim == 2 and idx not in total_mark]

gmsh.model.addPhysicalGroup(2, upper_mark, tag=5)
gmsh.model.addPhysicalGroup(2, lower_mark, tag=7)
gmsh.model.addPhysicalGroup(2, background_mark, tag=9)

upper_boundary = gmsh.model.getBoundary([(2, e) for e in upper_mark], recursive=False, oriented=False)
lower_boundary = gmsh.model.getBoundary([(2, e) for e in lower_mark], recursive=False, oriented=False)
outer_boundary = gmsh.model.getBoundary([(2, e) for e in background_mark], recursive=False, oriented=False)

upper_interface = [idx for (dim, idx) in upper_boundary if dim == 1]
lower_interface = [idx for (dim, idx) in lower_boundary if dim == 1]
total_interface = upper_interface + lower_interface
ext_boundary = [idx for (dim, idx) in outer_boundary if idx not in total_interface and dim == 1]

gmsh.model.addPhysicalGroup(1, upper_interface, tag=12)
gmsh.model.addPhysicalGroup(1, lower_interface, tag=15)
gmsh.model.addPhysicalGroup(1, ext_boundary, tag=17)

 # Create refined mesh around plates
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in upper_boundary+lower_boundary])

# Add threshold field for mesh refinement
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", plate_height/3)
gmsh.model.mesh.field.setNumber(2, "LcMax", plate_height*6)
gmsh.model.mesh.field.setNumber(2, "DistMin", plate_height*4)
gmsh.model.mesh.field.setNumber(2, "DistMax", plate_height*10)

gmsh.model.mesh.field.add("Min", 5)
gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
gmsh.model.mesh.field.setAsBackgroundMesh(5)

# Generate mesh
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
gmsh.write("capacitor.msh")

    
    

# Read mesh and set up connectivity
mesh_obj, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
mesh_obj.topology.create_connectivity(gdim - 1, gdim)  # Create necessary connectivity
gmsh.finalize()



# Save mesh for visualization
with XDMFFile(MPI.COMM_WORLD, "capacitor.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_obj)
    xdmf.write_meshtags(cell_tags, mesh_obj.geometry)


# Create function space
V = fem.functionspace(mesh_obj, ("CG", 2))


# Now locate DOFs based on the markers
upper_plate_dofs = fem.locate_dofs_topological(V, gdim - 1,facet_tags.find(12))
lower_plate_dofs = fem.locate_dofs_topological(V, gdim - 1, facet_tags.find(15))

# Apply Dirichlet boundary conditions
bc_upper = fem.dirichletbc(PETSc.ScalarType(V0), upper_plate_dofs, V)
bc_lower = fem.dirichletbc(PETSc.ScalarType(-V0), lower_plate_dofs, V)



# Combine boundary conditions
bcs = [bc_upper, bc_lower]    

# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Permittivity (εr = 1 for air)
epsilon_0 = 8.854e-12  # F/m
epsilon_r = 1.0        # relative permittivity for air
epsilon = epsilon_0 * epsilon_r

# Poisson equation: -∇·(ε∇u) = ρ/ε₀ (here ρ = 0 between plates)
a = epsilon * inner(grad(u), grad(v)) * dx
L = fem.Constant(mesh_obj, PETSc.ScalarType(0.0)) * v * dx

# Solve problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

# Save solution to XDMF format for visualization
with io.XDMFFile(mesh_obj.comm, "capacitor_solution.xdmf", "w") as file:
    file.write_mesh(mesh_obj)
    file.write_function(uh)


# Compute electric field (E = -grad(u))
gdim = mesh_obj.geometry.dim
E_field = -ufl.grad(uh)
V_E = fem.functionspace(mesh_obj, ("DG", 2, (gdim, )))
E = fem.Function(V_E)
E_expr = fem.Expression(E_field, V_E.element.interpolation_points())
E.interpolate(E_expr)


# Save the electric field
with io.VTXWriter(mesh_obj.comm, "capacitor_field.bp", E, engine="BP4") as vtx:
    vtx.write(0.0)

I am able to solve the plate capacitor problem. But the field lines on and near the plates are disoriented. I don’t it’s from paraview or fenics. Can someone please help?

Topic		Replies	Views
Misunderstanding regarding locate_dofs_topological General	2	334	March 1, 2023
A bug appears when using locate_dofs_topological(...) General boundary-conditions	1	81	April 22, 2024
Cannot find DOFS for discontinuous element dolfinx	8	470	April 12, 2023
Unable to properly locate dofs dolfinx	3	511	March 30, 2021
Using facet tags to define boundary conditions - GMSH - FEniCSx - dolfinx	4	1714	October 3, 2022

How to locate edges of a rectangle for defining dofs?

Related topics