...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
http://people.bath.ac.uk/em459/

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