1,848
edits
m (Remove redundant Category:Packages.) |
|||
(7 intermediate revisions by 3 users not shown) | |||
Line 8: | Line 8: | ||
== Tutorials == | == Tutorials == | ||
The solution of a problem can be logically divided in two steps. According to convenience or personal preference, they can be addressed with different files or just in one Octave script: | |||
* | * the description of the '''abstract problem''': this should be done via the Unified Form Language ('''UFL'''), which is ''a domain specific language for defining discrete variational forms and functionals in a notation '''close to pen-and-paper formulation'''.'' UFL is easy to learn, and in any case the User manual provides explanations and examples. [http://fenicsproject.org/documentation/ufl/1.2.0/user/user_manual.html#ufl-user-manual] As mentioned before, the problem can be defined in a separate .ufl file or handled directly in an m-file using ufl blocks. | ||
* a | * the implementation of a '''specific problem''', an instance of the abstract one: this is done in a script file ('''.m''') where the fem-fenics functions are used and the problem is solved. Their '''syntax is as close as possible to the python interface''', so that Fenics users should be comfortable with it, but it is also quite intuitive for beginners. The examples below show the equivalence between the different programming languages. | ||
Line 443: | Line 443: | ||
ufl J = derivative(F, u, du) | ufl J = derivative(F, u, du) | ||
ufl end | ufl end | ||
problem = 'HyperElasticity'; | problem = 'HyperElasticity'; | ||
Line 588: | Line 586: | ||
using namespace dolfin; | using namespace dolfin; | ||
Line 1,063: | 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,073: | 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. | ||
*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 '' | * 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,107: | 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 | [[Category:Octave Forge]] |