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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
"""This demo program solves Poisson's equation
- div C grad u(x, y) = f(x, y)
on the unit square with source f given by
f(x, y) = 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)
and boundary conditions given by
u(x, y) = 0 for x = 0 or x = 1
du/dn(x, y) = 0 for y = 0 or y = 1
The conductivity C is a symmetric 2 x 2 matrix which
varies throughout the domain. In the left part of the
domain, the conductivity is
C = ((1, 0.3), (0.3, 2))
and in the right part it is
C = ((3, 0.5), (0.5, 4))
The data files where these values are stored are generated
by the program generate_data.py
This demo is dedicated to BF and Marius... ;-)
"""
# Copyright (C) 2009-2011 Anders Logg
#
# 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: 2009-12-16
# Last changed: 2011-06-28
# Begin demo
from dolfin import *
# Read mesh from file and create function space
mesh = Mesh("../unitsquare_32_32.xml.gz")
V = FunctionSpace(mesh, "Lagrange", 1)
# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
# Define conductivity components as MeshFunctions
c00 = MeshFunction("double", mesh, "../unitsquare_32_32_c00.xml.gz")
c01 = MeshFunction("double", mesh, "../unitsquare_32_32_c01.xml.gz")
c11 = MeshFunction("double", mesh, "../unitsquare_32_32_c11.xml.gz")
# Code for C++ evaluation of conductivity
conductivity_code = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
#include <dolfin/mesh/MeshFunction.h>
class Conductivity : public dolfin::Expression
{
public:
// Create expression with 3 components
Conductivity() : dolfin::Expression(3) {}
// Function for evaluating expression on each cell
void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
{
const uint cell_index = cell.index;
values[0] = (*c00)[cell_index];
values[1] = (*c01)[cell_index];
values[2] = (*c11)[cell_index];
}
// The data stored in mesh functions
std::shared_ptr<dolfin::MeshFunction<double>> c00;
std::shared_ptr<dolfin::MeshFunction<double>> c01;
std::shared_ptr<dolfin::MeshFunction<double>> c11;
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<Conductivity, std::shared_ptr<Conductivity>, dolfin::Expression>
(m, "Conductivity")
.def(py::init<>())
.def_readwrite("c00", &Conductivity::c00)
.def_readwrite("c01", &Conductivity::c01)
.def_readwrite("c11", &Conductivity::c11);
}
"""
c = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(),
c00=c00, c01=c01, c11=c11, degree=0)
C = as_matrix(((c[0], c[1]), (c[1], c[2])))
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
a = inner(C*grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Save solution in VTK format
file = File("poisson.pvd")
file << u
# Plot solution
import matplotlib.pyplot as plt
plot(u)
plt.show()
|
