Fem-fenics: Difference between revisions

Jump to navigation Jump to search

Fem-fenics (edit)

Revision as of 19:47, 12 September 2014

4,959 bytes added ,  12 September 2014
→‎Tutorials: Add obstacle example
(→‎Relevant Implementation Details: Mention MPI parallelism)
(→‎Tutorials: Add obstacle example)
Line 1,059: Line 1,059:
# Hold plot
# Hold plot
interactive()
interactive()
© Copyright 2011, The FEniCS Project
</syntaxhighlight>}}
  </div>
</div>
<div style="clear:both"></div>
=== Obstacles in the Domain ===
[[File:Fem_fenics_Subdomains.png|right|500px]]
<math> - \mathrm {div} (a \nabla u) = 1 \quad \mbox { in } \Omega, </math>
<math> u = 5 \quad \mbox { on } \Gamma_T, </math>
<math> u = 0 \quad \mbox { on } \Gamma_B, </math>
<math> \nabla u \cdot \mathbf {n} = - 10 e^{- (y - 0.5)^2} \quad \mbox { on } \Gamma_L, </math>
<math> \nabla u \cdot \mathbf {n} = 1 \quad \mbox { on } \Gamma_R </math>
The above example is a weighted Poisson problem on the unit square. The diffusion coefficient <math> a </math> assumes value 0.01 on the obstacle <math> \Omega_1 = [0.5, 0.7] \times [0.2, 1.0] </math>, whilst 1 outside on <math> \Omega_0 = \Omega \setminus \Omega_1 </math>. You can find a detailed explanation on the [http://fenicsproject.org/documentation/dolfin/1.4.0/python/demo/documented/subdomains-poisson/python/documentation.html FEniCS website]. On the side, you can see the solution.
<div style="clear: both"></div>
<div style="width: 100%;">
  <div style="float:left; width: 48%">
{{Code|Define Poisson problem with obstacle with fem-fenics|<syntaxhighlight lang="octave" style="font-size:13px">
pkg load fem-fenics msh
# Create mesh
x = y = linspace (0, 1, 65);
[msh, facets] = Mesh (msh2m_structured_mesh (x, y, 0, 4:-1:1));
ufl start Subdomains
ufl fe = FiniteElement "(""CG"", triangle, 2)"
ufl u = TrialFunction (fe)
ufl v = TestFunction (fe)
ufl
ufl a0 = Coefficient (fe)
ufl a1 = Coefficient (fe)
ufl g_L = Coefficient (fe)
ufl g_R = Coefficient (fe)
ufl f = Coefficient (fe)
ufl
ufl a= "inner(a0*grad(u),grad(v))*dx(0) + inner(a1*grad(u),grad(v))*dx(1)"
ufl L = g_L*v*ds(1) + g_R*v*ds(3) + f*v*dx(0) + f*v*dx(1)
ufl end
V = FunctionSpace ("Subdomains", msh);
# Define problem coefficients
a0 = Constant ("a0", 1.0);
a1 = Constant ("a1", 0.01);
g_L = Expression ("g_L", @(x, y) - 10*exp(- (y - 0.5) ^ 2));
g_R = Constant ("g_R", 1.0);
f = Constant ("f", 1.0);
# Define subdomains
domains = MeshFunction ("dx", msh, 2, 0);
obstacle = SubDomain (@(x,y) (y >= 0.5) && (y <= 0.7) && ...
                      (x >= 0.2) && (x <= 1.0), false);
domains = mark (obstacle, domains, 1);
# Define boundary conditions
bc1 = DirichletBC (V, @(x, y) 5.0, facets, 2);
bc2 = DirichletBC (V, @(x, y) 0.0, facets, 4);
# Define variational form
a = BilinearForm ("Subdomains", V, V, a0, a1, domains);
L = LinearForm ("Subdomains", V, g_L, g_R, f, facets, domains);
# Solve problem
[A, b] = assemble_system (a, L, bc1, bc2);
sol = A \ b;
u = Function ("u", V, sol);
# Plot solution
[X, Y] = meshgrid (x, y);
U = u (X, Y);
surf (X, Y, U);
</syntaxhighlight>}}
  </div>
  <div style="float:right; width: 48%">
{{Code|Define Poisson problem with obstacle with FEniCS python |<syntaxhighlight lang="python" style="font-size:13px">
from dolfin import *
# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 0.0)
class Right(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 1.0)
class Bottom(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 0.0)
class Top(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 1.0)
class Obstacle(SubDomain):
  def inside(self, x, on_boundary):
    return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# Define mesh
mesh = UnitSquareMesh(64, 64)
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 2)
u = TrialFunction(V)
v = TestFunction(V)
# Define input data
a0 = Constant(1.0)
a1 = Constant(0.01)
g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))")
g_R = Constant("1.0")
f = Constant(1.0)
# Initialize mesh function for interior domains
obstacle = Obstacle()
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, 5.0, boundaries, 2),
      DirichletBC(V, 0.0, boundaries, 4)]
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]
# Define variational form
F = (inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)
    - g_L*v*ds(1) - g_R*v*ds(3)
    - f*v*dx(0) - f*v*dx(1))
# Separate left and right hand sides of equation
a, L = lhs(F), rhs(F)
# Solve problem
u = Function(V)
solve(a == L, u, bcs)
# Plot solution
plot(u, title="u")
interactive()
© Copyright 2011, The FEniCS Project
© Copyright 2011, The FEniCS Project
</syntaxhighlight>}}
</syntaxhighlight>}}
Eg123
22

edits

Retrieved from "https://wiki.octave.org/Special:MobileDiff/5051"

Navigation menu