Hello Community,
I am still relative new to FEniCS and this is my first forum post.
I was already able to create meshes with GMSH and import them into FEniCS using the msh to xdmf conversion method. Furthermore I was able to use the boundarys specified in GMSH to define DirechletBC in FEniCS. But I can’t figure out how to apply periodic BCs with the Line Tags created with GMSH.
My simple geo file:
// Mesh resulution
lc = 1e-1;
// Corner points of a rectangle
Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc};
Point(3) = {1, 1, 0, lc};
Point(4) = {0, 1, 0, lc};
// Lines in between the corner points
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
// Curve Loop and Surface of the rectangle
Curve Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
// Line Tags for BCs. Numeric names for FEniCS
Physical Curve(1111) = {4};
Physical Curve(4444) = {1};
Physical Curve(5555) = {3};
Physical Curve(9999) = {2};
Physical Surface(555) = {1};
Mesh conversion script:
import meshio
import numpy as np
# --------------------------------------------------
# ----------MESH CONVERSION----------
# Directory and name of the mesh which will be transformed
Meshdirectory = 'OrgMeshes/'
MeshFile = Meshdirectory + 'rect_periodic.msh'
def create_mesh(mesh, cell_type, prune_z=False):
# Read all relevant cells and get their physical data
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
# Read all points of the mesh
outpoints = mesh.points
if prune_z:
# In this case the z axis is deleted from the mesh (gmsh and salome always create 3d meshes)
outpoints = np.delete(arr=mesh.points, obj=2, axis=1)
# The relevant points and cells are combined to a meshio mesh and returned
return meshio.Mesh(points=outpoints, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]})
# Read in the mesh file
msh = meshio.read(MeshFile)
# Extract the mesh and the boundaries from the msh file
mesh_2d = create_mesh(msh, "triangle", True)
boundaries = create_mesh(msh, "line", True)
# Write the extracted mesh and boundaries to xdmf files (the h5 files are created automatically)
meshio.write("mesh.xdmf", mesh_2d)
meshio.write("boundaries.xdmf", boundaries)
# --------------------------------------------------
Mesh loading script:
import os
import csv
import matplotlib.pyplot as plt
import numpy as np
from dolfin import cpp
from fenics import *
# --------------------------------------------------
# ----------MESH LOADING----------
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
with XDMFFile("boundaries.xdmf") as infile:
infile.read(mvc, "name_to_read")
loadedBoundaryList = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=loadedBoundaryList, subdomain_id=12)
plot(mesh)
# --------------------------------------------------
# --------------------------------------------------
# ----------SIMULATION----------
D = 1 # Diffusion coefficient
T = 1.0 # Final time
num_steps = 100 # Number of time steps
num_steps_initial = 10 # Number of initial time steps
dt = T / num_steps # Time step size
dt_initial = dt / 10 # Initial time step size
dir_name = 'SimTest' # Folder for the results
V = FunctionSpace(mesh, 'P', 2)
INLET = DirichletBC(V, Constant(1), loadedBoundaryList, 1111)
OUTLET = DirichletBC(V, Constant(0), loadedBoundaryList, 9999)
bcs = [INLET, OUTLET]
# Define initial value
C_0 = interpolate(Constant(0), V)
# Define variational problem
C = TrialFunction(V)
v = TestFunction(V)
F_BE = (C - C_0) * v * dx + dt_initial * D * dot(grad(C), grad(v)) * dx
F_C = (C - C_0) * v * dx + 0.5 * dt * D * dot(grad(C), grad(v)) * dx + 0.5 * dt * D * dot(grad(C_0), grad(v)) * dx
F_BE_l, F_BE_r = lhs(F_BE), rhs(F_BE)
F_C_l, F_C_r = lhs(F_C), rhs(F_C)
How can i change the last script so that i not only apply the DirichletBCs for 1111 and 9999 but also a periodicBC between 4444 and 5555?
I hope I gave enough informations and my question is understandable otherwise please don’t hesitate to correct me or ask further questions.
Thank you all in advance,
Alex