helmholtz_resonator.py
3ed8cef19cfcmaster

helmholtz_resonator.py
View Options

			from matplotlib import pyplot as plt
			import fenics as fen
			import mshr
			import ufl
			from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE

			p = plt.jet()

			n = 50

			boundary = mshr.Polygon([fen.Point(0, 0), fen.Point(6, 0), fen.Point(6, 1),
			fen.Point(1.3, 1), fen.Point(1.3, 1.2),
			fen.Point(1.65, 1.2), fen.Point(1.65, 2.2),
			fen.Point(.65, 2.2), fen.Point(.65, 1.2),
			fen.Point(1, 1.2), fen.Point(1, 1), fen.Point(0, 1)])
			mesh = mshr.generate_mesh(boundary, n)

			for k in range(1):
			kappa = .25 + .05 * k
			ZR = 10.
			x, y = fen.SpatialCoordinate(mesh)[:]

			solver = HPE(kappa)
			solver.V = fen.FunctionSpace(mesh, "P", 3)
			solver.RobinBoundary = (lambda x, on_b: on_b and (fen.near(x[0], 6.)
			or fen.near(x[1], 2.2)))
			solver.NeumannBoundary = "REST"
			solver.signR = + 1.
			solver.NeumannDatum = [fen.Constant(0.),
			ufl.conditional(ufl.And(ufl.gt(y, 0.), ufl.lt(y, 1.)),
			fen.Constant(kappa), fen.Constant(0.))]
			solver.RobinDatumH = [fen.Constant(0.),
			ufl.conditional(ufl.gt(y, 1.25),
			fen.Constant(kappa / ZR),
			fen.Constant(kappa))]

			uh = solver.solve(kappa)[0]
			solver.plot(uh, name = "k={}".format(kappa))
			print(solver.norm(uh))

			from rrompy.hfengines.base import MarginalProxyEngine

			ms = MarginalProxyEngine(solver, [None])
			print(ms.A(1))