I don't see the boundary conditions in time

variational formulation
1

I have advection-diffusion problem, and solve test for square, i dot’h have my bc in time. My code:

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
plot(mesh)

dim = mesh.topology().dim()
mvc = MeshValueCollection("size_t", mesh, dim-1) 
with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
facet_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

F_dg = FunctionSpace(mesh, 'DG', 1)
V = VectorFunctionSpace(mesh, 'CG', 1)

bc_inlet = DirichletBC(F_dg, Constant(1), facet_domains, 7)
bc_w = DirichletBC(F_dg, Constant(0), facet_domains, 10)
bct = [bc_inlet, bc_w]
...
rig = lhs(F)
Lig = rhs(F)

gh = Function(F_dg)

A4 = assemble(rig)
[bc.apply(A4) for bc in bct]

gh.assign(t_n)

t = 0.0

for n in range(100):
    b4 = assemble(Lig)
    [bc.apply(b4) for bc in bct]
    
    solve(A4, gh.vector(), b4)
    plot(gh)
    
    max_t = np.abs(np.array(gh.vector())).max()
    print(max_t)
    t_n.assign(gh)

I can give full code if you need.

2

Please follow the guidelines: Read before posting: How do I get my question answered? making a minimal code example that reproduces the behaviour. Please read the second post in that thread as well to see good example of how a minimal working example should look like.

3

I thought the error was pretty simple, but for some reason I can’t see it, sorry
This minimal work code, where dont work boundary conditions:

from dolfin import * 

mesh = UnitSquareMesh(100, 100)
F_dg = FunctionSpace(mesh, 'DG', 1)
V = VectorFunctionSpace(mesh, 'CG', 1)

def boundary(x):
    return x[0] < 0.1 and x[1] < 0.1
bc = DirichletBC(F_dg, Constant(1), boundary)

u = Function(V)
u = project(Constant((1,1)), V)

t_n = Function(F_dg)
b = Expression("(x[0]>0.2 and x[0]<0.5) and (x[1] > 0.5 and x[1] < 0.8) ? 1.0 : 0.0", degree=1)
t_n = interpolate(b, F_dg)

plot(t_n)

t = TrialFunction(F_dg)
ti = TestFunction(F_dg)

Time = 1
num_times = 1000
dt = Time/num_times
k = Constant(dt)

alpha = Constant(2e1)
theta = Constant(1.0)
D = Constant(0.02)
f = Constant(0)

n = FacetNormal(mesh)
h = 2*Circumradius(mesh)

h_avg = (h('+') + h('-'))/2

un = (dot(u, n) + abs(dot(u, n)))/2.0

def a(x,z) :
    a_int = dot(grad(z), D*grad(x) - u*x)*dx

    a_fac = D('+')*(alpha('+')/h('+'))*dot(jump(z, n), jump(x, n))*dS \
        - D('+')*dot(avg(grad(z)), jump(x, n))*dS \
        - D('+')*dot(jump(z, n), avg(grad(x)))*dS


    a_vel = dot(jump(z), un('+')*x('+') - un('-')*x('-') )*dS \
        + dot(z, un*x)*ds

    a = a_int + a_fac + a_vel
    return a

a0=a(t_n,ti)
a1=a(t,ti)

A = (1/dt)*inner(t, ti)*dx - (1/dt)*inner(t_n,ti)*dx + theta*a1 + (1-theta)*a0

F = A

rig = lhs(F)
Lig = rhs(F)

gh = Function(F_dg)

A4 = assemble(rig)
bc.apply(A4)

gh.assign(t_n)

t = 0.0
for n in range(10):
    b4 = assemble(Lig)
    bc.apply(b4)    
    solve(A4, gh.vector(), b4)
    plot(gh)
    
    max_t = np.abs(np.array(gh.vector())).max()
    print(max_t)
    t_n.assign(gh)

I see my square but not see my boundaty conditions

4

This is because you are using a DG space. See: Dirichlet Boundary Conditions for Discontinuous Galerkin (DG) Methods - #4 by nate and the subsequent posts on how to apply Dirichlet conditions for such methods.

1 Like
5

Thanks, i have a quations, what’s mean “geometric” in:

bc = DirichletBC(V_dg, g, DirichletBoundary(), "geometric")

Thank you, I managed to solve the test problem, how can I set such conditions on the grids imported from gmsh?

6

There are usually two ways of locating degrees of freedom for setting boundary conditions:

  • geometrical (Check if a dof coordinate fulfills a certain condition)
  • topological (check if a facet fullfill a certain condition, then find all degrees of freedom related to the facet).

A big difference between CG and DG spaces is that DG degrees of freedom has no relation to global facets (as they only relate to a given cell).

To set such conditions for meshes from Gmsh, i suggest using a MeshFunction (imported through a mesh value collection) as described in for instance: Need help converting GMSH to FEniCS - #17 by dokken