...looking through Lawrence’s emails on the 3d solver, I found the code below, which seems to build the zhat field. It looks pretty impenetrable to me, is there any easier way of doing this or do I only need parts of it? Thanks, Eike ### BEGIN NORMALS BLOCK ### # Build new extruded coordinate function space normal_fs = VectorFunctionSpace(mesh, "DG", 0) zhat = Function(normal_fs) base_coords = mesh._old_mesh.coordinates base_coords_dim = base_coords.function_space().cdim base_coords_arity = base_coords.cell_node_map().arity expr = Expression(('x[0]', 'x[1]', 'x[2]')) v0 = lambda x: ast.Symbol("v0", (x,)) v1 = lambda x: ast.Symbol("v1", (x,)) n = lambda x: ast.Symbol("n", (x,)) x = lambda x: ast.Symbol("x", (x,)) coords = lambda x, y: ast.Symbol("base_coords", (x, y)) normals = lambda x, y: ast.Symbol("normals", (x, y)) body = [] body += [ast.Decl("double", v(3)) for v in [v0, v1, n, x]] body.append(ast.Decl("double", "dot")) body.append(ast.Decl("double", "norm")) body.append(ast.Assign("norm", 0.0)) body.append(ast.Assign("dot", 0.0)) body.append(ast.Decl("int", "i")) body.append(ast.Decl("int", "c")) body.append(ast.Decl("int", "d")) body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), [ast.Assign(v0("i"), ast.Sub(coords(1, "i"), coords(0, "i"))), ast.Assign(v1("i"), ast.Sub(coords(2, "i"), coords(0, "i"))), ast.Assign(x("i"), 0.0)])) # n = v0 x v1 body.append(ast.Assign(n(0), ast.Sub(ast.Prod(v0(1), v1(2)), ast.Prod(v0(2), v1(1))))) body.append(ast.Assign(n(1), ast.Sub(ast.Prod(v0(2), v1(0)), ast.Prod(v0(0), v1(2))))) body.append(ast.Assign(n(2), ast.Sub(ast.Prod(v0(0), v1(1)), ast.Prod(v0(1), v1(0))))) body.append(ast.For(ast.Assign("i", 0), ast.Less("i", 3), ast.Incr("i", 1), [ast.Incr(x(j), coords("i", j)) for j in range(3)])) body.extend([ast.FlatBlock("dot += (%(x)s) * n[%(i)d];\n" % {"x": x, "i": i}) for i, x in enumerate(expr.code)]) body.extend([ast.FlatBlock("norm += n[%(i)d] * n[%(i)d];\n" % {"i": i}) for i in range(3)]) body.append(ast.FlatBlock("norm = sqrt(norm);\n")) # normals[d][c] = n[c] / norm; assign_block = [ast.Assign(normals("d", "c"), ast.Div(n("c"), ast.Symbol("norm")))] body.append(ast.IMul("norm", ast.Ternary(ast.Less("dot", 0), "-1", "1"))) body.append(ast.For(ast.Assign("d", 0), ast.Less("d", 1), ast.Incr("d", 1), [ast.For(ast.Assign("c", 0), ast.Less("c", base_coords_dim), ast.Incr("c", 1), assign_block)])) # Looks complicated, actually just produces the code below # void normal_calc (double** base_coords , double** normals ) # { # double v0[3] ; # double v1[3] ; # double n[3] ; # double x[3] ; # double dot ; # double norm ; # norm = 0.0; # dot = 0.0; # int i ; # int c ; # int d ; # for (i = 0; i < 3; ++i) # { # v0[i] = base_coords[1][i] - base_coords[0][i]; # v1[i] = base_coords[2][i] - base_coords[0][i]; # x[i] = 0.0; # } # n[0] = v0[1] * v1[2] - v0[2] * v1[1]; # n[1] = v0[2] * v1[0] - v0[0] * v1[2]; # n[2] = v0[0] * v1[1] - v0[1] * v1[0]; # for (i = 0; i < 3; ++i) # { # x[0] += base_coords[i][0]; # x[1] += base_coords[i][1]; # x[2] += base_coords[i][2]; # } # dot += (x[0]) * n[0]; # dot += (x[1]) * n[1]; # dot += (x[2]) * n[2]; # norm += n[0] * n[0]; # norm += n[1] * n[1]; # norm += n[2] * n[2]; # norm = sqrt(norm); # norm *= (dot < 0 ? -1 : 1); # for (d = 0; d < 1; ++d) # { # for (c = 0; c < 3; ++c) # { # normals[d][c] = n[c] / norm; # } # } # } kernel = op2.Kernel(ast.FunDecl("void", "normal_calc", [ast.Decl("double**", "base_coords"), ast.Decl("double**", "normals")], ast.Block(body)), "normal_calc") base_coords = mesh._old_mesh.coordinates op2.par_loop(kernel, zhat.cell_set, base_coords.dat(op2.READ, base_coords.cell_node_map()), zhat.dat(op2.WRITE, zhat.cell_node_map())) ### END NORMALS BLOCK ###

On 5 Dec 2014, at 10:01, Eike Mueller <E.Mueller@bath.ac.uk> wrote:

Dear firedrakers,

I have the following function spaces:

W_2 = HDiv(U_2 x V_0) + HDiv(U_1 x V_1) [velocity] W_2^b = U_2 x V_0 [buoyancy]

where U_0 -> U_1 -> U_2 and V_0 -> V_1 are horizontal- and vertical de Rham complexes.

I need to be able to apply the operator Q: W_2^b -> W_2 defined by

Q_{ij} = <\vec{w}_i,\hat{z} \gamma_j>

where \vec{w}_i and \gamma_j are basis functions and \hat{z} is the unit normal in the vertical direction. To apply Q to a buoyancy field b, I want to write something like

u = assemble(dot(TestFunction(W_2),zhat*b)*dx)

but how do I get the zhat in there? Is that even possible?

Thanks a lot,

Eike

--

Dr Eike Hermann Mueller Research Associate (PostDoc)

Department of Mathematical Sciences University of Bath Bath BA2 7AY, United Kingdom

+44 1225 38 5803 e.mueller@bath.ac.uk <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/

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