Errors in use of conditional expressions in variational form

variational formulation
10

Thank you very much. This seems to have resolved my issue. I also had to create the function p using conditional as opposed to expression.

#p_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
#                   degree = 3)

x0, x1 = MeshCoordinates(mesh)

p_Ex = conditional(ge(0, x0-x1), -2*(exp(x0-x1)-1)/(exp(x0-x1)-1), -(x0-x1))

My functioning code in full now reads as follows:

from fenics import *
from dolfin import *
from mshr import *
import numpy as np
#import ufl                                                                     #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/


#ufl.algorithms.apply_derivatives.CONDITIONAL_WORKAROUND = True                 #from https://fenicsproject.org/qa/10763/how-can-create-an-expression-is-dependent-on-function-value/    

mesh = UnitSquareMesh(16, 16)

#***************Defining variables and corresponding function spaces**********

W = FiniteElement('DG', mesh.ufl_cell(), 0)                                    #Trial/test space of Discontinuous Lagrange elements  
V = FiniteElement('RT', mesh.ufl_cell(), 1)                                    #Trial/test space of Raviart-Thomas elements
element = MixedElement([W,V])
M = FunctionSpace(mesh, element)                                               #Creating mixed finite element space
W_0 = M.sub(0).collapse()                                                      #Retrieving trial/test space of constant elements 
V_0 = M.sub(1).collapse()                                                      #Retrieving trial/test space of vector elements
    
    
#**************************boundary condition*********************************     

#p_Ex = Expression('x[0]-x[1] >= 0 ? -2*(exp(x[0]-x[1])-1)/(exp(x[0]-x[1])-1): -(x[0]-x[1])',
#                   degree = 3)

x0, x1 = MeshCoordinates(mesh)

p_Ex = conditional(ge(0, x0-x1), -2*(exp(x0-x1)-1)/(exp(x0-x1)-1), -(x0-x1))

def boundary(x, on_boundary):
    return on_boundary

#bc = DirichletBC(M.sub(0), u_Ex, boundary)

#bc = DirichletBC(W_0, u_Ex, boundary)


#*******************Initial Condition******************************************

p = project(p_Ex, W_0)                                                         #Projection of u_Ex onto scalar component of mixed space.

bc = DirichletBC(M.sub(0), p, boundary)

#*****************Setting values of certain parameters************************

k = 2.0                                                                        # Linearisation constant

T = 0.14                                                                       #Final time in corresponding parabolic problem
N = 10                                                                         #Number of steps in parabolic problem
tau = T/N                                                                      #Time step size in parabolic problem


#*****************Defining variational form************************************

#uq = Function(M)                                                               

uq = TrialFunction(M)                                                          #Defining trial functions of the mixed formulation

#u, q = uq.split(True)

u,q = split(uq)

w,v = TestFunctions(M)                                                         #Defining test functions of the mixed formulation

bn_ = conditional(le(0, p), pow(pi,2)/2 - pow(p,2)/2, pow(pi,2)/2)             #The function b acting on the projection u_p of u_Ex 

a1 = k*(u*w*dx) + tau*((div(q))*w*dx)                                          #bilinear form in equation for u 

L1 = (bn_ + k*p)*w*dx                                                          #linear form in equation for u 
  
a2 = dot(q, v)*dx - u*(div(v))*dx                                              #bilinear form in equation for q

L2 = 0                                                                         #linear form in equation for q  

a = a1 + a2

L = L1 + L2     

#*****************Solving mixed problem****************************************

uq = Function(M)                                                               #Solution

solve(a==L, uq, bc)                                                            #Solving mixed variational form 
        
u, q = uq.split(True)                                                          #Splitting solution into functions from the respective components of the mixed space 

print(u)

unodes = u.compute_vertex_values(mesh)                                         #Storing vertex values of numerical solution

print(unodes)                                                                  #Displaying vertex values of numerical solution

Thanks again Nate and David :slight_smile: .

2 Likes