Re: [firedrake] Tetrahedral mesh
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_...
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
Hi Andrew, The mesh turned out as I needed. However, solving the 3d Helmholtz equation on the same mesh could not converge to a solution. Could you inspect the following code:
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 <mailto: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 <mailto: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 ofhttps://github.com/firedrakeproject/firedrake/blob/master/firedrake/utility_... <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 <mailto: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 <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 10/12/15 14:12, Angwenyi David wrote:
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
This is the indefinite form of the Helmholtz equation. Did you actually want the sign-positive Helmholtz equation? The later has strong form: -laplacian u + u = f When integrating by parts the weak form becomes: <grad u, grad v> + <u, v> = <f, v> As to incorporation of boundary conditions, see http://firedrakeproject.org/variational-problems.html#incorporating-boundary... Cheers, Lawrence
Integration by parts is equivalent to using Green’s function; and so it turns out as follows: laplacian u = - u + f Greens function: integ((laplacian u)v) + integ(dot(grad u, grad v)) = surfinteg((partial_n u)v) So that the weak form is obtained by substituting laplacian u = -u + f from the Helmholtz equation in the Green’s function. The same is arrived at by integration by parts. Sincerely, David
On 10 Dec 2015, at 15:23, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
On 10/12/15 14:12, Angwenyi David wrote:
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
This is the indefinite form of the Helmholtz equation. Did you actually want the sign-positive Helmholtz equation?
The later has strong form:
-laplacian u + u = f
When integrating by parts the weak form becomes:
<grad u, grad v> + <u, v> = <f, v>
As to incorporation of boundary conditions, see
http://firedrakeproject.org/variational-problems.html#incorporating-boundary...
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
Andrew McRae
-
Angwenyi David
-
Lawrence Mitchell