Solving for power dissipation where k*grad(u)**2.0 How to solve for this?

Question

Hello All. I am new to Fenics and have been having some issues. I am trying to solve for power dissipation that occurs due to joule heating. I start by solving standard poisson problem with 2 different materials, and then attempt to solve for k*grad(u)**2.0 after I solve for u. However this leads to some odd results where my power dissipation can becomes negative at the material boundary change. This is physically impossible ( and should not happen in the solution either as the grad(u) term is squared. Where am I going wrong here?

I attached an example of the code `

from dolfin import *

import sys, math, numpy

nx = 32; ny = 32
mesh = UnitSquareMesh(nx, ny)

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

class Omega0(SubDomain):
def inside(self, x, on_boundary):
return True if x[1] <= 0.5 else False

class Omega1(SubDomain):
def inside(self, x, on_boundary):
return True if x[1] >= 0.5 else False

subdomain0 = Omega0()
subdomain0.mark(subdomains, 0)
subdomain1 = Omega1()
subdomain1.mark(subdomains, 1)

V0 = FunctionSpace(mesh, 'DG', 0)
k = Function(V0)

print 'mesh:', mesh
print 'subdomains:', subdomains
print 'k:', k

k_values = [1.5, 50] # values of k in the two subdomains
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
k.vector()[cell_no] = k_values[subdomain_no]

help = numpy.asarray(subdomains.array(), dtype=numpy.int32)
k.vector()[:] = numpy.choose(help, k_values)

print 'k degree of freedoms:', k.vector().array()

V = FunctionSpace(mesh, 'Lagrange', 1)

tol = 1E-14 # tolerance for coordinate comparisons
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1]) < tol

Gamma_0 = DirichletBC(V, Constant(0), BottomBoundary())

class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] - 1) < tol

Gamma_1 = DirichletBC(V, Constant(1), TopBoundary())

bcs = [Gamma_0, Gamma_1]

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = kinner(nabla_grad(u), nabla_grad(v))dx
L = fvdx
u = Function(V)
solve(a == L, u, bcs)

disipation=project(k*grad(u)**2.0,FunctionSpace(mesh,'Lagrange',1))

file = File("disipation.pvd")
file << disipation

Comments

I did not run your code, but I guess you observe some artifacts, which may be largest at the material jump, right? Did you try to project your result to DG0 instead? This is actually the space where k*grad(u)**2 makes sense. Projecting piecewise constants to H1 functions can lead to weird effects.

---

Thank you. That does fix the problem. Switching to second order lagrange function also removed the negative power dissipation issue, but results in a smooth transition between the two materials which is expected but is accurate on low resolution mesh. I need to read up on the spaces in Fenics and how they work.

Thank you again!