60
edits
Fem-fenics: Difference between revisions
Jump to navigation
Jump to search
no edit summary
No edit summary |
No edit summary |
||
| Line 118: | Line 118: | ||
<div style="clear:both"></div> | <div style="clear:both"></div> | ||
=== Mixed Formulation for Poisson Equation === | |||
<div style="width: 100%;"> | |||
<div style="float:left; width: 48%"> | |||
{{Code|Define Poisson problem|<syntaxhighlight lang="octave" style="font-size:13px"> | |||
pkg load fem-fenics msh | |||
import_ufl_Problem ('MixedPoisson') | |||
# Create mesh | |||
x = y = linspace (0, 1, 33); | |||
mesh = Mesh(msh2m_structured_mesh (x, y, 1, 1:4)); | |||
# File MixedPoisson.ufl | |||
# BDM = FiniteElement("BDM", triangle, 1) | |||
# DG = FiniteElement("DG", triangle, 0) | |||
# W = BDM * DG | |||
V = FunctionSpace('MixedPoisson', mesh); | |||
# Define trial and test function | |||
# File MixedPoisson.ufl | |||
# (sigma, u) = TrialFunctions(W) | |||
# (tau, v) = TestFunctions(W) | |||
# CG = FiniteElement("CG", triangle, 1) | |||
# f = Coefficient(CG) | |||
f = Expression ('f', @(x,y) 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)); | |||
# Define variational form | |||
# File MixedPoisson.ufl | |||
# a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx | |||
# L = - f*v*dx | |||
a = BilinearForm ('MixedPoisson', V); | |||
L = LinearForm ('MixedPoisson', V, f); | |||
# Define essential boundary | |||
bc1 = DirichletBC (SubSpace (V, 1), @(x,y) [0; -sin(5.0*x)], 1); | |||
bc2 = DirichletBC (SubSpace (V, 1), @(x,y) [0; sin(5.0*x)], 3); | |||
# Compute solution | |||
[A, b] = assemble_system (a, L, bc1, bc2); | |||
sol = A \ b; | |||
func = Function ('func', V, sol); | |||
sigma = Function ('sigma', func, 1); | |||
u = Function ('u', func, 2); | |||
# Plot solution | |||
plot (sigma); | |||
plot (u); | |||
</syntaxhighlight>}} | |||
</div> | |||
<div style="float:right; width: 48%"> | |||
{{Code|Define Poisson problem|<syntaxhighlight lang="python" style="font-size:13px"> | |||
from dolfin import * | |||
# Create mesh | |||
mesh = UnitSquareMesh(32, 32) | |||
# Define function spaces and mixed (product) space | |||
BDM = FunctionSpace(mesh, "BDM", 1) | |||
DG = FunctionSpace(mesh, "DG", 0) | |||
W = BDM * DG | |||
# Define trial and test functions | |||
(sigma, u) = TrialFunctions(W) | |||
(tau, v) = TestFunctions(W) | |||
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") | |||
# Define variational form | |||
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx | |||
L = - f*v*dx | |||
# Define function G such that G \cdot n = g | |||
class BoundarySource(Expression): | |||
def __init__(self, mesh): | |||
self.mesh = mesh | |||
def eval_cell(self, values, x, ufc_cell): | |||
cell = Cell(self.mesh, ufc_cell.index) | |||
n = cell.normal(ufc_cell.local_facet) | |||
g = sin(5*x[0]) | |||
values[0] = g*n[0] | |||
values[1] = g*n[1] | |||
def value_shape(self): | |||
return (2,) | |||
G = BoundarySource(mesh) | |||
# Define essential boundary | |||
def boundary(x): | |||
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS | |||
bc = DirichletBC(W.sub(0), G, boundary) | |||
# Compute solution | |||
w = Function(W) | |||
solve(a == L, w, bc) | |||
(sigma, u) = w.split() | |||
# Plot sigma and u | |||
plot(sigma) | |||
plot(u) | |||
interactive() | |||
</syntaxhighlight>}} | |||
</div> | |||
</div> | |||
<div style="clear:both"></div> | |||
[[Category:OctaveForge]] | [[Category:OctaveForge]] | ||
[[Category:Packages]] | [[Category:Packages]] | ||