|
Hi All, I'm just getting started with Dolfin/FFC and have modifed several of the demos into a working python code for a test problem, but am having difficulty getting the equivalent cpp code working. The variational form looks like ---------------------------------------------------------------- element = FiniteElement("Lagrange", "triangle", 2) v = TestFunction(element) u = TrialFunction(element) f = Function(element) # set domain size coefficient h/delta hondel = 64. c = 1/hondel/hondel a = f*f*f*dot(grad(v),grad(u))*dx + c*f*v*u*dx L = v.dx(1)*f*f*f*dx ------------------------------------------------------------------ (which arises in a non-linear wave problem) where f is a user defined function (porosity) which needs to get passed to the bilinear form a. The python code (attached) works just fine, but in the c++ code where I create ------------------------------------------------- // Set up problem UnitSquare mesh(32, 32); // Create source function Source f(mesh); // Dirichlet boundary conditions Function u1(mesh,1.0); DirichletBoundary boundary; DirichletBC bc(u1, mesh, boundary); // Solution function Function u; // Create forms and nonlinear PDE PressureEquationBilinearForm a(u); PressureEquationLinearForm L(f); LinearPDE pde(a, L, mesh, bc); // Solve pde.solve(u); ---------------------------------------------------------------- Everthing compiles just fine but at runtime I get... Initializing DOLFIN version 0.7.3. Initializing PETSc with given command-line arguments. Creating linear PDE with one boundary condition. Solving linear PDE. terminate called after throwing an instance of 'std::runtime_error' what(): *** Error: Function contains no data. I'm sure I'm just being stupid, but all help greatly appreciated cheers marc |
# This demo program calculates the pressureEquation potential # in the unit square by solving # # - div f^3 grad u + c fu = -div f^3 \khat # # where f is the porosity and c is a constant related to the box height in compaction lengths # # The boundary conditions are: # # u(x, y) = 1 for y = 0, y=1 # du(x,y) = 0 for x=0,x=1
__author__ = "Marc Spiegelman ([EMAIL PROTECTED])"
__date__ = "2008-02-14"
__copyright__ = ""
__license__ = "GNU LGPL Version 2.1"
from dolfin import *
hondel = 64; # size of domain in compaction lengths
c = hondel*hondel
# Create mesh and finite element
mesh = UnitSquare(64,64)
# Finite elements
element = FiniteElement("Lagrange", "triangle", 2) # last argument: order of the polynomial
# Porosity term
class Porosity(Function):
def __init__(self, element, mesh):
Function.__init__(self, element, mesh)
def eval(self, values, x):
dx = x[0] - 0.5
dy = x[1] - 0.5
sigma = 0.1
values[0] = 1.+3.*exp(-(dx*dx + dy*dy)/sigma/sigma)
# Dirichlet boundary condition
class DirichletFunction(Function):
def __init__(self, element,mesh):
Function.__init__(self, element,mesh)
def eval(self, values, x):
values[0] = 0.0
# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((abs(x[1] - 1.) < DOLFIN_EPS or abs(x[1]) < DOLFIN_EPS) and on_boundary and True)
# Define variational problem
v = TestFunction(element)
u = TrialFunction(element)
f = Porosity(element,mesh)
a = f*f*f*dot(grad(v),grad(u))*dx + c*f*v*u*dx
L = v.dx(1)*f*f*f*dx
# Define boundary condition
u0 = DirichletFunction(element,mesh)
boundary = DirichletBoundary()
bc = DirichletBC(u0, mesh, boundary)
# Solve PDE and plot solution
pde = LinearPDE(a, L, mesh, bc)
u = pde.solve()
plot(u)
plot(f)
# Save solution to file
f1 = File("pressure.pvd")
f1 << u
f2 = File("porosity.pvd")
f2 << f
#Hold plot
interactive()
PressureEquation.form
Description: Binary data
main.cpp
Description: Binary data
---------------------------------------------------- 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
