-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbackend2.py
152 lines (122 loc) · 5.87 KB
/
backend2.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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
from mpi4py import MPI
import dolfinx
import dolfin
import pyvista
from dolfinx.fem import (VectorFunctionSpace,
FunctionSpace,
Function,
Constant,
locate_dofs_topological,
dirichletbc, Expression,
assemble_scalar, form)
from dolfinx.mesh import locate_entities_boundary, compute_midpoints
from ufl import (TrialFunction,
TestFunction, Measure, dot,
sym, grad, dx, inner, tr,
nabla_div, Identity, sqrt)
from mshr import Polygon, generate_mesh
from dolfinx.fem.petsc import LinearProblem
import numpy as np
class Simulation:
def __init__(self, corners, meshDensity = 15) -> None:
self.corners = corners
self.domain = self.createPolygonalMesh(corners, meshDensity)
self.V = VectorFunctionSpace(self.domain, ("Lagrange", 1))
self.U = FunctionSpace(self.domain, ("DG", 0))
self.bcs = []
def createPolygonalMesh(self, corners, meshDensity=15):
corners = list(map(dolfin.Point, corners))
domain = Polygon(corners)
domain = generate_mesh(domain, meshDensity)
with dolfin.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf") as file:
file.write(domain)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh()
self.locs = compute_midpoints(domain, 2, np.arange(domain.topology.index_map(2).size_local, dtype='int32'))
return domain
def createFunctions(self):
self.u = TrialFunction(self.V)
self.v = TestFunction(self.V)
self.B = Function(self.V, name="Body Force")
self.T = Constant(self.domain, dolfinx.default_scalar_type((0, 0)))
self.density = Function(self.U, name="Density")
self.vm = Function(self.U, name="Von Mises Stress")
self.complianceArr = Function(self.U, name="Compliance")
def fixedJoint(self, f):
fdim = self.domain.topology.dim - 1
boundary_facets = locate_entities_boundary(self.domain, fdim, f)
u_D = np.array([0, 0], dtype=dolfinx.default_scalar_type)
dofs = locate_dofs_topological(self.V, fdim, boundary_facets)
print(len(dofs))
self.bcs.append(dirichletbc(u_D, dofs, self.V))
def rollingJoint(self, f):
fdim = self.domain.topology.dim - 1
boundary_facets = locate_entities_boundary(self.domain, fdim, f)
dofs = locate_dofs_topological(self.V.sub(1), fdim, boundary_facets)
self.bcs.append(dirichletbc(dolfinx.default_scalar_type(0), dofs, self.V.sub(1)))
def applyForce(self, forces):
def bodyForce(x):
T = np.zeros_like(x[:2,:])
for loc, force in forces:
c = np.isclose(x[0], loc[0], rtol=0.01) & np.isclose(x[1], loc[1], rtol=0.01)
T = T + np.stack([force[0]*c, force[1]*c], axis=0)
return T
self.B.interpolate(bodyForce)
def constituentEqns(self, Emin, E0, penal, nu):
E = Emin + self.density**penal*(E0-Emin)
lambda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
self.eps = lambda u: sym(grad(u))
self.sigma = lambda u: lambda_ * nabla_div(u) * Identity(len(u)) + 2 * mu * self.eps(u)
s = lambda u: self.sigma(u) - 1. / 3 * tr(self.sigma(u)) * Identity(len(u))
self.von_Mises = lambda u: sqrt(3. / 2 * inner(s(u), s(u)))
def solve(self):
ds = Measure('ds', domain=self.domain)
a = inner(self.sigma(self.u), self.eps(self.v)) * dx
L = dot(self.B, self.v) * dx + dot(self.T, self.v) * ds
problem = LinearProblem(a, L, bcs=self.bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
self.disp = problem.solve()
def updateStress(self):
stress_expr = Expression(self.von_Mises(self.disp), self.U.element.interpolation_points())
comp_expr = Expression(inner(self.eps(self.disp), self.sigma(self.disp)), self.U.element.interpolation_points())
self.complianceArr.interpolate(comp_expr)
self.vm.interpolate(stress_expr)
def compliance(self):
comp = MPI.COMM_WORLD.allreduce(assemble_scalar(form(inner(self.eps(self.disp), self.sigma(self.disp)) * dx)),op=MPI.SUM)
return comp
def displayResult(self, showUndeformed = False):
if pyvista.OFF_SCREEN:
pyvista.start_xvfb(wait=0.1)
p = pyvista.Plotter()
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(self.domain, self.domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
if showUndeformed:
p.add_mesh(grid)
grid["u"] = self.disp.x.array.reshape((geometry.shape[0], 2))
z = np.zeros(len(grid['u']))
grid['u'] = np.c_[grid["u"], z]
warped = grid.warp_by_vector("u", factor=1)
warped.cell_data["VonMises"] = self.vm.vector.array
warped.set_active_scalars("VonMises")
p.add_mesh(warped,clim=[0, 0.3*max(self.vm.vector.array)], cmap='jet')
p.add_camera_orientation_widget()
if not pyvista.OFF_SCREEN:
p.show()
p.screenshot('imgs/stress.png')
else:
print("Unable to show plot.")
def displayDistribution(self):
if pyvista.OFF_SCREEN:
pyvista.start_xvfb(wait=0.1)
p = pyvista.Plotter()
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(self.domain, self.domain.topology.dim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
grid.cell_data["Density"] = 1-self.density.vector.array
grid.set_active_scalars("Density")
p.add_mesh(grid, cmap='gray', clim=[0, 1])
p.add_camera_orientation_widget()
if not pyvista.OFF_SCREEN:
p.show()
p.screenshot('imgs/density.png')
else:
print("Unable to show plot.")