This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

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

0 votes

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

closed with the note: The question was answered by Christian Waluga, use DG0 space instead to get proper result.
asked Apr 26, 2016 by soviet911 FEniCS Novice (180 points)
closed Apr 27, 2016 by soviet911

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!

...