dolfin team mailing list archive

[mspieg <mspieg@xxxxxxxxxxxxxxxxx>]

"""This demo attempts to project velocity and pressure gradients from a previously calculated Stoke's problem
pressure onto a conservative fluid-velocity space

i.e. Given discrete functions V (in P2) and p (P1) from a Stokes problem.
Form a new vector valued function

    v = V - grad(p)

by projecting V and grad(p) onto the function space describing v.  Note, all three discrete
functions share the same mesh 

Modified from the PyDolfin mixed-poisson demo.py by Kristian B Oelgaard
"""

__author__    = "Marc Spiegelman"
__date__      = " 2 Mar 2009 17:31:54"
__copyright__ = "Copyright (C) 2009 Marc Spiegelman"
__license__   = "GNU LGPL Version 2.1"

from dolfin import *

# read in Velocity and pressure functions (from a Taylor-Hood Stokes problem)
p = Function("pressure.xml.gz")
V = Function("velocity.xml.gz")

# extract the mesh
mesh = p.function_space().mesh()
#mesh = V.function_space().mesh()

# define new function space for v on the same mesh
BDM = FunctionSpace(mesh, "BDM", 1)

# Define projection problem
w = TestFunction(BDM)
u = TrialFunction(BDM)

a = dot(w,u)*dx

# various linear forms (some work, some don't)

# project both V and grad(p)...this is desired but broken
L =  dot(w,V)*dx - dot(w,grad(p))*dx

# project just grad(p), this works if mesh is extracted from p
#L =  -dot(w,grad(p))*dx 

# project just V. this works if mesh is extracted from V
# L =  dot(w,V)*dx # works if mesh is extracted from V


# Compute solution
problem = VariationalProblem(a, L)
v = problem.solve()

# Project v for post-processing
P1 = VectorFunctionSpace(mesh, "CG", 1)
v_proj = project(v, P1)

# Plot solution
plot(v_proj)
interactive()


# Save solution to file
f0 = File("melt_velocity.xml")
f0 << v

# Save solution to pvd format
f1 = File("melt_velocity.pvd")
f1 << v_proj