-
Notifications
You must be signed in to change notification settings - Fork 1
/
main00_alternative_mesh.py
121 lines (102 loc) · 5.4 KB
/
main00_alternative_mesh.py
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
# -*- coding: utf-8; -*-
"""Generate a very simplified uniform mesh for testing, and tag its boundaries."""
import typing
import numpy as np
import matplotlib.pyplot as plt
import dolfin
from extrafeathers import autoboundary
from extrafeathers import meshiowrapper
from extrafeathers import meshmagic, plotmagic
from .config import mesh_filename, Boundaries, L, N, aspect
def main():
assert dolfin.MPI.comm_world.size == 1, "Mesh can only be generated in serial mode, please run without mpirun."
# https://fenicsproject.org/olddocs/dolfin/latest/cpp/de/dac/classdolfin_1_1UnitSquareMesh.html#a0cff260683cd6632fa569ddda6529a0d
# 'std::string diagonal ("left", "right", "right/left", "left/right", or "crossed")
# indicates the direction of the diagonals.'
# mesh = dolfin.UnitSquareMesh(N, N)
# mesh = dolfin.UnitSquareMesh(N, N, "crossed")
# mesh = dolfin.UnitSquareMesh.create(N, N, dolfin.CellType.Type.quadrilateral)
# mesh = meshmagic.trimesh(N, N) # rows of equilateral triangles
# mesh = meshmagic.trimesh(N, N, align="y") # columns of equilateral triangles
# from dolfin import ALE, Constant
# ALE.move(mesh, Constant((-0.5, -0.5)))
# Aspect ratio (= width / height) of domain.
# `N` should be integer-divisible by this.
def vtxpreproc(vtxs):
def grade(xx, dense_side, strength): # xx: np.array of points in [0, 1]
assert dense_side in ("left", "right")
# # root, singular at the coarse side
# if dense_side == "right":
# return xx**(1 / strength)
# return 1.0 - (1.0 - xx)**(1 / strength) # left
# # exponential, smooth everywhere
if dense_side == "right":
return 1.0 - (np.exp(strength * xx) - 1.0) / (np.exp(strength * 1) - 1.0)
return (np.exp(strength * xx) - 1.0) / (np.exp(strength * 1) - 1.0) # left
def mean(a, b):
return a + (b - a) / 2
# # Grade the mesh. In the 3D printing test case, most stuff happens at the upper left corner.
# # With load ramping, this is optional.
# vtxs[:, 0] = grade(vtxs[:, 0], "left", 2)
# vtxs[:, 1] = grade(vtxs[:, 1], "right", 2)
# # # vtxs[:, 0] = mean(vtxs[:, 0], grade(vtxs[:, 0], "left", 2))
# # # vtxs[:, 1] = mean(vtxs[:, 1], grade(vtxs[:, 1], "right", 2))
# # vtxs[:, 0] = mean(grade(vtxs[:, 0], "left", 2), grade(vtxs[:, 0], "left", 2))
# # vtxs[:, 1] = mean(grade(vtxs[:, 1], "right", 2), grade(vtxs[:, 1], "right", 2))
# center mesh on origin
vtxs[:, 0] -= 0.5
vtxs[:, 1] -= 0.5
# scale uniformly to desired domain length
vtxs *= L
# scale `y` to account for desired aspect ratio
vtxs[:, 1] /= aspect
return vtxs
mesh = meshmagic.trimesh(nx=N, ny=N // aspect, align="y", vtxpreproc=vtxpreproc)
# mesh = meshmagic.trimesh(nx=N, ny=N // aspect, align="x", vtxpreproc=vtxpreproc)
domain_parts = dolfin.MeshFunction('size_t', mesh, mesh.topology().dim(), mesh.domains())
if mesh.cell_name() == "quadrilateral":
V = dolfin.FunctionSpace(mesh, "Q", 1)
else:
V = dolfin.FunctionSpace(mesh, "P", 1)
ignored_cells, nodes_dict = meshmagic.all_cells(V)
ignored_dofs, nodes_array = meshmagic.nodes_to_array(nodes_dict)
xmin = np.min(nodes_array[:, 0])
xmax = np.max(nodes_array[:, 0])
ymin = np.min(nodes_array[:, 1])
ymax = np.max(nodes_array[:, 1])
print(f"Generated mesh with {nodes_array.shape[0]} vertices, {len(ignored_cells)} cells.")
# `autoboundary` calls the callback for each facet on the external boundary
# (i.e. the facet is on a boundary, and no neighboring subdomain exists).
def autoboundary_callback(submesh_facet: dolfin.Facet, fullmesh_facet: dolfin.Facet) -> typing.Optional[int]:
p = submesh_facet.midpoint()
x, y = p.x(), p.y()
on_vert_boundary = dolfin.near(y, ymin) or dolfin.near(y, ymax)
on_horz_boundary = dolfin.near(x, xmin) or dolfin.near(x, xmax)
if dolfin.near(x, xmin) and not on_vert_boundary:
return Boundaries.LEFT.value
elif dolfin.near(x, xmax) and not on_vert_boundary:
return Boundaries.RIGHT.value
elif dolfin.near(y, ymax) and not on_horz_boundary:
return Boundaries.TOP.value
elif dolfin.near(y, ymin) and not on_horz_boundary:
return Boundaries.BOTTOM.value
return None # this facet is not on a boundary we are interested in
# Tag the boundaries.
print("Tagging boundaries...")
boundary_parts: dolfin.MeshFunction = autoboundary.find_subdomain_boundaries(fullmesh=mesh,
submesh=mesh,
subdomains=domain_parts,
boundary_spec={},
callback=autoboundary_callback)
print("Writing HDF5 file...")
meshiowrapper.write_hdf5_mesh(mesh_filename, mesh, None, boundary_parts)
print("Showing mesh.")
plotmagic.mpiplot_mesh(V)
plt.axis("equal")
plotmagic.plot_facet_meshfunction(boundary_parts, names=Boundaries)
plt.legend(loc="best")
plt.title("Generated mesh")
plt.show()
print("All done.")
if __name__ == "__main__":
main()