>>> nx = 10

>>> ny = 10

>>> nz = 20

>>> mesh = BoxMesh(nx, ny, nz, 1.0, 1.0, 2.0)

>>> V = FunctionSpace(mesh, "CG", 1)

>>> u = TrialFunction(V)

>>> v = TestFunction(V)

>>> f = Function(V)

>>> f.interpolate(Expression("(2*x[0])+(3*x[1])+(0.5*x[2])"))

Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', tetrahedron, 1, dim=3), 3), FiniteElement('Lagrange', tetrahedron, 1)), 4)

>>> a = (dot(grad(v), grad(u)) - v * u) * dx

>>> L = -f * v * dx

>>> u = Function(V)

>>> solve(a == L, u, solver_parameters={'ksp_type': 'cg’})

>>> File("helmholtz.pvd") << u

>>> output = File("helmholtz.pvd")
>>> output<<u
>>> del output

I will be glad also to get a comment on how to incorporate boundary terms.

Thanks in advance.

Sincerely,

David

On 09 Dec 2015, at 16:23, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:

This is after doing mesh = BoxMesh (...)

On 9 Dec 2015 15:19, "Angwenyi David" <kipkoej@gmail.com> wrote:

Hi Andrew,

Here is what I get:

>>> File("mesh.pvd")<<mesh.coordinates

Traceback (most recent call last):

File "<stdin>", line 1, in <module>

AttributeError: 'module' object has no attribute ‘coordinates'

What could be the issue?

On 09 Dec 2015, at 15:36, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:

Hi David,

The mesh for this is BoxMesh. For documentation, type '? BoxMesh' at a python prompt, after importing firedrake, or see line 475 onwards of https://github.com/firedrakeproject/firedrake/blob/master/firedrake/utility_meshes.py

In your case, you want

BoxMesh(nx, ny, nz, 1.0, 1.0, 2.0), where nx, ny and nz are the number of cells in the x, y and z direction respectively. Ideally (see below), each cuboid is divided into 6 tetrahedra, so the total number of cells is 6*nx*ny*nz, so you can control this exactly.

Unfortunately, there is a trap: if you want cells whose x, y and z dimensions are very different, the mesh generator will generate extra points.

With the dimensions you gave above, if you have nx = 10, ny = 10, nz = 20, you get what you expect - lots of cubes, each divided into 6 tetrahedra.

However, if you have nx = 5, ny = 5, nz = 100, you get something completely different (lots of small cells).

You can visualise the mesh by using

File("mesh.pvd") << mesh.coordinates

and opening the file in Paraview in 'Wireframe' mode.

There's a very easy workaround, so do reply if this is a problem for you.

Best,

Andrew

On 9 December 2015 at 14:15, Angwenyi David <kipkoej@gmail.com> wrote:

Dear all,

I would like to solve a 3d Helmholtz equation in a domain 0 <= x <= 1, 0 <= y <= 1, 0 <= z <= 2 tessellated into any number of tetrahedra. Is this possible within Firedrake, and if yes, what is the mesh name?

David
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake