(Deleted) Can facet jump() be used on a Function()?


68
views
0
7 months ago by
I have a time stepping algorithm where I need to use the jump on an interface labeled 1 from a previous solve for u0 to update formulate the jump condition on the current solve. This is from implementing a simple ODE on the boundary of the form, like  $g_t=\Delta t\left[\left[u_{t-1}\right]\right]^2+g_{t-1}$gt=Δt[[ut1]]2+gt1  in the example. However, in the code below, jump(u0) appears to be equivalent to u0 when the arrays are compared. Is jump( ) able to be used on a function, such as u0, that is the result of a previous DG solution time step?
from fenics import *
import numpy as np

ofilename = 'model.xml'
mesh = Mesh(ofilename)
subdomains = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

V = FunctionSpace(mesh, 'DG', 2)
u = TrialFunction(V)
v = TestFunction(V)

n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2
alpha = Constant(10.0)
gamma = Constant(alpha)

p0 = Constant(1.0)
g0 = interpolate(Constant(1.0), V)
g1 = Function(V)
g2 = Function(V)

time = np.linspace(0.0, 1.0, 3.0)
for t in range(len(time)-1):
    dt = time[t+1] - time[t]
    F_interior = inner(grad(v), grad(u))*dx
    F_edges = -inner(jump(v,n), avg(grad(u)))*dS \
              - inner(avg(grad(v)), jump(u,n))*dS \
              + alpha/h_avg*inner(jump(u), jump(v))*dS  
              
    F_bc = -inner(grad(v), (u-p0)*n)*ds(4) \
            - inner(v*n, grad(u))*ds(4) \
            + (gamma/h)*v*(u-p0)*ds(4)
            
    F_bc += -inner(grad(v), (u+p0)*n)*ds(6) \
             - inner(v*n, grad(u))*ds(6) \
             + (gamma/h)*v*(u+p0)*ds(6)   
             
    F = F_interior + F_edges + F_bc       
    a = lhs(F)
    l = rhs(F)
    u0 = Function(V)
    solve(a==l, u0)
    
    g1.interpolate(project(dt*jump(u0)*jump(u0) + g0, V))
    g2.interpolate(project(dt*u0*u0 + g0, V))
    print(g1.vector().array()[:])
    print(g2.vector().array()[:])
    
    g0.assign(g1)
Community: FEniCS Project
Please login to add an answer/comment or follow this question.
The thread is closed. No new answer/comment may be added.

Similar posts:
Search »
  • Nothing matches yet.