dolfin team mailing list archive

[mspieg <mspieg@xxxxxxxxxxxxxxxxx>]

# 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 (mspieg@xxxxxxxxxxxxxxxxx)"
__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()

Attachment: PressureEquation.form
Description: Binary data

Attachment: main.cpp
Description: Binary data