|
Hi all, I've recently been having some issues when trying to solve with forms where I attach multiple discrete functions that share the same mesh but not necessarily the same function_space (i.e. element). This used to work under 0.8.1 and I suspect it's a bug. A reasonably simple example of the problem can be stated as follows. let V, p be discrete Velocity and pressure functions generated from a Stoke's solve with Taylor-hood mixed elements. The two solutions share the same mesh, but are in different function spaces (V is P2, p is P1). I now want to form a new vector valued function v = V - grad(p) by projecting V and grad(p) onto an appropriate function space. (If you're wondering, this is a toy version of a deformable porous media problem) Taking what seems to be the obvious approach ---------------------------------------------- # 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() # define a 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 L = dot(w,V)*dx - dot(w,grad(p))*dx # Compute solution problem = VariationalProblem(a, L) v = problem.solve() ----------------------------------------------------------- generates the following error from solve (JIT/FFC seems to work fine) Solving linear variational problem Traceback (most recent call last): File "test.py", line 44, in <module> v = problem.solve() File "/Users/mspieg/mspieg/hg_repositories/FEniCS/srcFEniCS/dolfin/local/lib/python2.5/site-packages/dolfin/variationalproblem.py", line 52, in solve cpp.VariationalProblem.solve(self, u) File "/Users/mspieg/mspieg/hg_repositories/FEniCS/srcFEniCS/dolfin/local/lib/python2.5/site-packages/dolfin/cpp.py", line 10742, in solve return _cpp.VariationalProblem_solve(*args) RuntimeError: *** Error: Unable to extract mesh from form (nonmatching meshes for function spaces). If I just project the pressure part of v i.e. L = -dot(w,grad(p)) everything works fine. Similarly if I just project the velocity part (and extract the mesh from V rather than p) i.e. mesh = V.function_space().mesh() ... L = dot(w,V) that also works fine. If i use the mesh from p to project V, things break and vice-versa (and yet it should be the same mesh. And they certainly are if you write them out and diff the .xml files) I can obviously project each part separately and add them together, but that seems like a hack and I don't see why this shouldn't work. I also see this behaviour in my c++ codes as well and have tracked it to problems where I try to attach multiple discrete functions on different function spaces as coefficients to a form. Like I said, this used to work in 0.8.1, but breaks now. In general, I don't see why you can't have arbitrary discrete functions in a linear form as long as they can be evaluated during assembly. apologies for the long e-mail but all help greatly appreciated cheers marc p.s. I've attached a test script with the above problem and variations p.p.s versions should all be up to date development versions Dolfin 0.9.1 dev ufc 1.1.1 dev FFC 0.6.1 dev instant 0.9.6 dev (changset 283) |
"""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
---------------------------------------------------- Marc Spiegelman Lamont-Doherty Earth Observatory Dept. of Applied Physics/Applied Math Columbia University tel: 845 704 2323 (SkypeIn) ---------------------------------------------------- |
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
