helmholtz_resonator.py
No OneTemporary
Actions

Subscribers

None

File Metadata

	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))