1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
|
"""This demo demonstrates how to solve a mixed Poisson type equation
defined over a sphere (the surface of a ball in 3D) including how to
create a cell_orientation map, needed for some forms defined over
manifolds."""
# Copyright (C) 2012 Marie E. Rognes
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# First added: 2012-12-09
# Last changed: 2012-12-09
# Begin demo
from dolfin import *
import numpy
import matplotlib.pyplot as plt
# Read mesh
mesh = Mesh("../sphere_16.xml.gz")
# Define global normal
global_normal = Expression(("x[0]", "x[1]", "x[2]"), degree=1)
mesh.init_cell_orientations(global_normal)
# Define function spaces and basis functions
RT1 = FiniteElement("RT", mesh.ufl_cell(), 1)
DG0 = FiniteElement("DG", mesh.ufl_cell(), 0)
R = FiniteElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement((RT1, DG0, R)))
(sigma, u, r) = TrialFunctions(W)
(tau, v, t) = TestFunctions(W)
g = Expression("sin(0.5*pi*x[2])", degree=2)
# Define forms
a = (inner(sigma, tau) + div(sigma)*v + div(tau)*u + r*v + t*u)*dx
L = g*v*dx
# Tune some factorization options
if has_petsc():
# Avoid factors memory exhaustion due to excessive pivoting
PETScOptions.set("mat_mumps_icntl_14", 40.0)
PETScOptions.set("mat_mumps_icntl_7", "0")
# Avoid zero pivots on 64-bit SuperLU_dist
PETScOptions.set("mat_superlu_dist_colperm", "MMD_ATA")
# Solve problem
w = Function(W)
solve(a == L, w, solver_parameters={"symmetric": True})
(sigma, u, r) = w.split()
# Plot CG1 representation of solutions
sigma_cg = project(sigma, VectorFunctionSpace(mesh, "CG", 1))
u_cg = project(u, FunctionSpace(mesh, "CG", 1))
plt.figure()
plot(sigma_cg)
plt.figure()
plot(u_cg)
plt.show()
# Store solutions
file = File("sigma.pvd")
file << sigma
file = File("u.pvd")
file << u
|
