Fem-fenics: Difference between revisions

Jump to navigation Jump to search

Fem-fenics (edit)

Revision as of 11:12, 10 June 2019

5,295 bytes added ,  10 June 2019
m
Remove redundant Category:Packages.
(→‎Tutorials: Mention ufl blocks in the opening clauses)
m (Remove redundant Category:Packages.)
 
(5 intermediate revisions by 3 users not shown)
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 *
# Define mesh
mesh = UnitSquareMesh(64, 64)
# Create classes for defining parts of the boundaries
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)
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
# 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
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)))
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>}}
Line 1,069: Line 1,303:
The relevant implementation details which the user should know are:
The relevant implementation details which the user should know are:


*all the objects are managed using boost::shared_ptr <>. It means that '''the same resource can be shared by more objects''' and useless copies should be avoided. For example, if we have two different functional spaces in the same problem, like with Navier-Stokes for the velocity and the pressure, the mesh is shared between them and no one has its own copy.
* All the objects are managed using boost::shared_ptr <>. It means that '''the same resource can be shared by more objects''' and useless copies should be avoided. For example, if we have two different functional spaces in the same problem, like with Navier-Stokes for the velocity and the pressure, the mesh is shared between them and no one has its own copy.


*The '''essential BC are imposed directly to the matrix''' with the command '''assemble()''', which sets to zero all the off diagonal elements in the corresponding line, sets to 1 the diagonal element and sets to the exact value the rhs. This means that we could loose the symmetry of the matrix, if any. To avoid this problem and preserve the symmetry of the system it is available the '''assemble_system()''' command which builds at once the lhs and the rhs.
* The '''essential BC are imposed directly to the matrix''' with the command '''assemble()''', which sets to zero all the off diagonal elements in the corresponding line, sets to 1 the diagonal element and sets to the exact value the rhs. This means that we could loose the symmetry of the matrix, if any. To avoid this problem and preserve the symmetry of the system it is available the '''assemble_system()''' command which builds at once the lhs and the rhs.


*The '''coefficient of the variational problem''' can be specified using either an ''Expression()'', a ''Constant()'' or a ''Function()''. They are different objects which behave in different ways: an ''Expession'' or a ''Constant'' object overloads the eval() method of the  dolfin::Expression class and it is evaluated at run time using the octave function feval (). A ''Function'' instead doesn't need to be evaluated because it is assembled copying element-by-element the values contained in the input vector.
* The '''coefficient of the variational problem''' can be specified using either an ''Expression()'', a ''Constant()'' or a ''Function()''. They are different objects which behave in different ways: an ''Expression'' or a ''Constant'' object overloads the eval() method of the  dolfin::Expression class and it is evaluated at run time using the octave function feval (). A ''Function'' instead doesn't need to be evaluated because it is assembled copying element-by-element the values contained in the input vector.


* Unfortunately the feature used in previous versions of the package to handle ''DirichletBC'' is not implemented yet in the FEniCS library for distributed execution. This means that prior to running code via MPI it should be adapted to use the newly introduced ''MeshFunction'' to mark border subsets, so that essential boundary conditions are correctly applied.


== Known issues ==
== Known issues ==
Line 1,103: Line 1,338:
* Presentation at MOX [https://drive.google.com/file/d/0ByWLfuWVSWHbZWZzRzY2em5PU28/edit?usp=sharing]
* Presentation at MOX [https://drive.google.com/file/d/0ByWLfuWVSWHbZWZzRzY2em5PU28/edit?usp=sharing]


[[Category:Octave-Forge]]
[[Category:Octave Forge]]
Retrieved from "https://wiki.octave.org/Special:MobileDiff/5047...11972"

Navigation menu