dolfin team mailing list archive

[Melanie Jahny <question155482@xxxxxxxxxxxxxxxxxxxxx>]

Question #155482 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/155482

    Status: Answered => Open

Melanie Jahny is still having a problem:
Thanks so much for your work. Now I can apply the robin boundary condition.
I still have one question left concerning the dirichlet boundary condition. When applying cell_domains = subdomains in A=assemble(a,...) I always got a matrix error:
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Matrix is missing diagonal entry in row 12860!

Do I have to fill the diagonal matrix entry with small values (10^-4 for example)?
How can I do that in FEniCS? Or do you think my matrix is wrong?

And is it right to set:
interior_facet_domains=boundary_parts, exterior_facet_domains=boundary_parts
in the "assemble" function? (sorry, I just want to be sure that my code is right)

In your first answer you asked me to give you my complete code, so here
it is:

#####################################################
### IMPLEMENTATION OF ADVEKTION DIFFUSION EQUATION
#####################################################


# --- DEFINE PARAMETERS -------------------------------------------
x_max = 154.0
y_max = 78.0 
membrane_max_y = 49.0 

# mesh parameters
x_div = 154 
y_div = 78

# --- define rectangle mesh ---------------------------------------
mesh = Rectangle(0,0,x_max,y_max,x_div, y_div)
mesh.order()

# time parameters
dt = 0.1
t  = 0.1
T  = 1.0

# Parameters for initial condition
v_mem = 0.002 # velocity at the membrane
Delta = 0.001


# Diffusion constants
D_const = 4.07*1e-11 # mm^2/s   

 # Function Space and Functions
V = FunctionSpace(mesh, "Lagrange", 1)
conc = TrialFunction(V) # solution
conc_prev = Function(V) # solution from previous time step
eta = TestFunction(V)

# Set initial condition
                    
c_init = Expression("exp(-%g*%g*(x[1]-%g)*(x[1]-%g)/(%g*%g))*exp(-0.5*((x[0]-0.1*%g)/(%g*0.001))*((x[0]-0.1*%g)/(%g*0.001)))" %\
                          (v_mem,v_mem, membrane_max_y, membrane_max_y, D_const, D_const, x_max, x_max, x_max, x_max ))

conc_prev.interpolate(c_init)

# ----------------------------------------------------------------------------------
velocity = Constant((0.0,-1.0))
# ----------------------------------------------------------------------------------

subdomains = MeshFunction('uint', mesh, 2) # 2 dimensional subdomains

# define channel where advection diffusion equation should be solved
class Channel(SubDomain):  
   def inside(self,x,on_boundary):
       return True if x[1] >= 0.0 else False
       
sub_channel = Channel()
sub_channel.mark(subdomains,0)


# --- --- linear form for time derivative ----------------------------------------
F_time = inner((conc-conc_prev),eta)*dx(0) 
   
# --- --- linear form for diffusion
F_Diff = D_const*inner(grad(conc), grad(eta))*dx(0)
   
# --- --- linear form for advection ---------------------------------------------   
F_Adv = inner(inner(velocity, grad(conc)),eta)*dx(0)


########################################################
# ------ set robin boundary condition
normal = FacetNormal(mesh)

boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

# robin boundary condition at membrane
class RobinBoundary_Mem(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] == membrane_max_y and on_boundary
         

Gamma_Mem = RobinBoundary_Mem()
Gamma_Mem.mark(boundary_parts,0)

class RobinBoundary_Left(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]== 0.0 and x[1] > membrane_max_y

# robin boundary condition at left side
Gamma_Left = RobinBoundary_Left()
Gamma_Left.mark(boundary_parts,1)

bound_mem = ((conc('+')*velocity('+')[0]*normal('+')[0] + conc('-')*velocity('-')[0]*normal('-')[0]) + \
                 (conc('+')*velocity('+')[1]*normal('+')[1] + conc('-')*velocity('-')[1]*normal('-')[1]))\
                   *eta('+')*dS(0) 
                   
bound_left = -(inner(conc*velocity,normal)*eta)*ds(1)

#---------------------------------------------------------------------------------------------

F = dt*(F_Adv + F_Diff + bound_mem + bound_left) + F_time

a = lhs(F) 
L = rhs(F)


#------------------------------------------------------------------------------------------
# --- --- Dirichlet Boundary Condition

c_const = Constant(0.0)

boundary_up = compile_subdomains('x[1] == 78.0 && on_boundary')
bc_up = DirichletBC(V, c_const, boundary_up)

boundary_right = compile_subdomains('x[0] == 154.0 && x[1] > 49.0 && on_boundary')   
bc_right = DirichletBC(V, c_const, boundary_right)    

bc = [bc_up, bc_right]

#-------------------------------------------------------------------------------------------------------

conc = conc_prev

# !!!!!!!!!!!!!!!!!
#A = assemble(a, cell_domains=subdomains, interior_facet_domains=boundary_parts, exterior_facet_domains=boundary_parts) # doesnt work
A = assemble(a, interior_facet_domains = boundary_parts, exterior_facet_domains=boundary_parts)
  
# time loop  
while(t<T):
     
   b = assemble(L, cell_domains = subdomains, interior_facet_domains=boundary_parts, exterior_facet_domains=boundary_parts) 
 
   bc_up.apply(A, b)
   bc_right.apply(A, b)      

   solve(A, conc.vector(),b)
   conc_prev.assign(conc)
      
   t += dt

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.