August 2016 - firedrake - imperial.ac.uk

Re: [firedrake] Boundary IDs on extruded meshes
by Andrew McRae 02 Aug '16

02 Aug '16

Hi Hanie, An ExtrudedMesh inherits boundary ids from the 'base mesh', plus the extra identifiers "bottom" and "top". So, the sides x==0 through y==1 are id 1-4, because this is what the base mesh has. Then z==0 is "bottom" and z==1 is "top". However, your code is failing because (I believe) you're trying to set one component of a VectorFunctionSpace to 0. Firedrake allows you to do this on normal meshes, but doesn't yet support this on extruded meshes, hence the error message. (You could set all components of the vector to 0 by using just wSpace.sub(0), but presumably this isn't what you want). There are various workarounds you could use for now, such as setting the boundary condition weakly, but they're not completely satisfying. Andrew On 1 August 2016 at 21:15, Hanie Joodat <hanie.joodat88(a)gmail.com> wrote: > Dear all, > > I want to enforce DirichletBCs on an extruded mesh. For a normal unit cube > tetrahedron mesh I could use the following face ID's: > > 1: x ==0 > 2: x==1 > 3: y==0 > 4: y==1 > 5: z==0 > 6: z==1 > > However, when I do the following: > > meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) > mesh = ExtrudedMesh(meshbase, Nz, Dz) > > I get this error: > > Traceback (most recent call last): > File "VMS_patch_test_3D_Double_firedrake.py", line 51, in <module> > FrontBackBC1 = DirichletBC(wSpace.sub(0).sub(1),Constant(0.0),(3,4)) > File "/data/Software/firedrake/src/firedrake/firedrake/bcs.py", line 68, > in __init__ > raise NotImplementedError("Indexed VFS bcs not implemented on extruded > meshes") > NotImplementedError: Indexed VFS bcs not implemented on extruded meshes > > What should I do to properly identify my boundary faces? > > Regards, > > Hanie >

2 1

Boundary IDs on extruded meshes
by Hanie Joodat 01 Aug '16

01 Aug '16

Dear all, I want to enforce DirichletBCs on an extruded mesh. For a normal unit cube tetrahedron mesh I could use the following face ID's: 1: x ==0 2: x==1 3: y==0 4: y==1 5: z==0 6: z==1 However, when I do the following: meshbase = RectangleMesh(Nx, Ny, Lx, Ly, quadrilateral=True) mesh = ExtrudedMesh(meshbase, Nz, Dz) I get this error: Traceback (most recent call last): File "VMS_patch_test_3D_Double_firedrake.py", line 51, in <module> FrontBackBC1 = DirichletBC(wSpace.sub(0).sub(1),Constant(0.0),(3,4)) File "/data/Software/firedrake/src/firedrake/firedrake/bcs.py", line 68, in __init__ raise NotImplementedError("Indexed VFS bcs not implemented on extruded meshes") NotImplementedError: Indexed VFS bcs not implemented on extruded meshes What should I do to properly identify my boundary faces? Regards, Hanie

1 0

Re: [firedrake] Unexpected asymmetry of solution of Poisson equation
by Andrew McRae 01 Aug '16

01 Aug '16

Miklos is exactly right. What you see in the very corners is the following: if you have a quad with 3 zeroes and a nonzero (let's say a 1), if the quad is split along a 0-0 diagonal the you get a 0-0-0 triangle and a 0-1-0 triangle, but if the quad is split a 0-1 diagonal then you get two 0-0-1 triangles. It's even worse if the quad is +1 -1 -1 +1 in which case you get a 'valley' or a 'ridge', even though the actual function is symmetric about 0. To answer your last question, the triangular mesh isn't symmetric in that way. On 1 August 2016 at 11:08, George Ovchinnikov <lives9(a)gmail.com> wrote: > Dear all, > > I'm solving Poisson equation and I found some irregularity. Here the code: > > from firedrake import * > > mesh = UnitSquareMesh(8, 8, quadrilateral=True) > > V = FunctionSpace(mesh, "CG", 1) > W = FunctionSpace(mesh, "DG", 0) > w = Function(W) > w.interpolate(Expression(1)) #Initial material distribution > > bc0 = DirichletBC(V, Expression(("0")), 1) > bc1 = DirichletBC(V, Expression(("0")), 2) > bc2 = DirichletBC(V, Expression(("0")), 3) > bc3 = DirichletBC(V, Expression(("0")), 4) > bcs = [bc0, bc1, bc2, bc3] > > u = TrialFunction(V) > v = TestFunction(V) > > f = Function(V) > f.interpolate(Expression("1")) > L = f * v * dx > rhs = assemble(L) > > a = (w * dot(grad(u), grad(v))) * dx > frm = assemble(a, bcs = bcs) > sol = Function(V) > solve(frm, sol, rhs, solver_parameters = {'ksp_type': 'preonly', > 'pc_type': 'lu'}) > File("u.pvd") << sol > > The asymmetry of solution (see the figure attached) is clearly visible: > for example, upper left element is not the same as upper right. > > I'm wondering what causes this and can this be remedied? > > Also, I'm getting exactly same result (with the same asymmetry) on > triangular, not quadrilateral mesh. Why solutions of triangular and > quadrilateral meshes are the same? > > Sincerely, > George >

2 1

Unexpected asymmetry of solution of Poisson equation
by George Ovchinnikov 01 Aug '16

01 Aug '16

Dear all, I'm solving Poisson equation and I found some irregularity. Here the code: from firedrake import * mesh = UnitSquareMesh(8, 8, quadrilateral=True) V = FunctionSpace(mesh, "CG", 1) W = FunctionSpace(mesh, "DG", 0) w = Function(W) w.interpolate(Expression(1)) #Initial material distribution bc0 = DirichletBC(V, Expression(("0")), 1) bc1 = DirichletBC(V, Expression(("0")), 2) bc2 = DirichletBC(V, Expression(("0")), 3) bc3 = DirichletBC(V, Expression(("0")), 4) bcs = [bc0, bc1, bc2, bc3] u = TrialFunction(V) v = TestFunction(V) f = Function(V) f.interpolate(Expression("1")) L = f * v * dx rhs = assemble(L) a = (w * dot(grad(u), grad(v))) * dx frm = assemble(a, bcs = bcs) sol = Function(V) solve(frm, sol, rhs, solver_parameters = {'ksp_type': 'preonly', 'pc_type': 'lu'}) File("u.pvd") << sol The asymmetry of solution (see the figure attached) is clearly visible: for example, upper left element is not the same as upper right. I'm wondering what causes this and can this be remedied? Also, I'm getting exactly same result (with the same asymmetry) on triangular, not quadrilateral mesh. Why solutions of triangular and quadrilateral meshes are the same? Sincerely, George

2 1

Re: [firedrake] Product in an expression
by David Ham 01 Aug '16

01 Aug '16

Hi Floriane, It looks rather like you are trying to express the basis function expansion at the Firedrake level. This is something which is likely to be quite hard because you're operating at the wrong abstraction level. Can you post at a weak form (i.e. PDEs and function spaces) level what it is you are trying to achieve? Many thanks, David On Mon, 1 Aug 2016 at 16:29 Floriane Gidel [RPG] <mmfg(a)leeds.ac.uk> wrote: > Dear all, > > > I would like to compute the following sum over i : phi(x,y,z,t) = > psi_i(x,y,t)varphi_i(z), where : > > - each psi_i (x,y,t) is a Firedrake function defined on a horizontal plane, > > -varphi_i is defined by a lagrange expansion of order n, that is: > varphi_i(z) = product( (z-z_k)/(z_i - z_k) ) for k = 0:n , and k different > from i. Each z_k can be computed analytically and denotes the discretised z > coordinate. For instance, z_k can be defined as z_k = (n-k)H0/n , for z in > [0,H0]. > > > How can I define varphi_i so that I can compute the dot product > dot(psi,varphi) on Firedrake ? I tried two options: > > - use the symbpy package to compute the expression of varphi_i with the > 'product' function. But in that case the expression of varphi_i depends on > z as a symbol, and I did not manage to use it correctly on Firedrake > afterwards, > > - define varphi_i(z) as a Firedrake function in z, but in that case, how > can I compute the product ? I would normally use an expression, but is > there a 'product expression' that I could interpolate ? Something like : > varphi[i].interpolate(Expression( "product((x[2]-z[k]) / (z[i]-z[k]) )", > z=z, i=i, for k=0:n)) > > > Thank you for your help, > > Floriane >

2 1

Product in an expression
by Floriane Gidel [RPG] 01 Aug '16

01 Aug '16

Dear all, I would like to compute the following sum over i : phi(x,y,z,t) = psi_i(x,y,t)varphi_i(z), where : - each psi_i (x,y,t) is a Firedrake function defined on a horizontal plane, -varphi_i is defined by a lagrange expansion of order n, that is: varphi_i(z) = product( (z-z_k)/(z_i - z_k) ) for k = 0:n , and k different from i. Each z_k can be computed analytically and denotes the discretised z coordinate. For instance, z_k can be defined as z_k = (n-k)H0/n , for z in [0,H0]. How can I define varphi_i so that I can compute the dot product dot(psi,varphi) on Firedrake ? I tried two options: - use the symbpy package to compute the expression of varphi_i with the 'product' function. But in that case the expression of varphi_i depends on z as a symbol, and I did not manage to use it correctly on Firedrake afterwards, - define varphi_i(z) as a Firedrake function in z, but in that case, how can I compute the product ? I would normally use an expression, but is there a 'product expression' that I could interpolate ? Something like : varphi[i].interpolate(Expression( "product((x[2]-z[k]) / (z[i]-z[k]) )", z=z, i=i, for k=0:n)) Thank you for your help, Floriane

1 0