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
|
"""This demo program demonstrates how to manipulate (higher-order) mesh
coordinates."""
# Copyright (C) 2016 Jan Blechta
#
# 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/>.
from dolfin import *
import matplotlib.pyplot as plt
# Create mesh
comm = MPI.comm_world
mesh = UnitDiscMesh.create(comm, 20, 2, 2)
plt.figure()
plot(mesh)
# Fetch coordinate function
C = FunctionSpace(mesh, mesh.ufl_coordinate_element())
c = Function(C)
get_coordinates(c, mesh.geometry())
# Deform coordinates harmonically subject to BC
u, v = TrialFunction(C), TestFunction(C)
a = inner(grad(u), grad(v))*dx
L = dot(Constant((0, 0)), v)*dx
bc1 = DirichletBC(C, (-1, -1), "x[0] < -0.5")
bc2 = DirichletBC(C, c, "x[0] >= -0.5")
displacement = Function(C)
solve(a == L, displacement, [bc1, bc2])
c_vec = c.vector()
c_vec += displacement.vector()
# Set coordinates
set_coordinates(mesh.geometry(), c)
plt.figure()
plot(mesh)
# We can create (cubic) mesh from function
C3 = VectorFunctionSpace(mesh, "Lagrange", 4)
c3 = interpolate(c, C3)
mesh3 = create_mesh(c3)
plt.figure()
plot(mesh3)
# Display plots
plt.show()
|
