1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
|
# Copyright (c) 2005 Johan Hoffman
# Licensed under the GNU GPL Version 2
#
# Modified by Anders Logg 2006
#
# First added: 2005
# Last changed: 2006-03-28
#
# The contiuity equation for the incompressible
# Navier-Stokes equations using cG(1)cG(1)
#
# Compile this form with FFC: ffc NSEContinuity2D.form.
name = "NSEContinuity3D"
scalar = FiniteElement("Lagrange", "tetrahedron", 1)
vector = VectorElement("Lagrange", "tetrahedron", 1)
constant_scalar = FiniteElement("Discontinuous Lagrange", "tetrahedron", 0)
cell = "tetrahedron"
K1 = VectorElement("Lagrange", cell, 1)
# Dimension of domain
d = K1.cell_dimension()
q = TestFunction(scalar) # test basis function
P = TrialFunction(scalar) # trial basis function
uc = Function(vector) # linearized velocity
delta1 = Function(constant_scalar) # stabilization parameter
P0 = Function(scalar) # trial basis function
k = Constant("tetrahedron") # time step
i0 = Index() # index for tensor notation
i1 = Index() # index for tensor notation
def ugradu(u, v):
return [dot(u, grad(v[i])) for i in range(d)]
alpha1 = 2.0 * k
# Bilinear and linear forms
a = alpha1*dot(grad(q), grad(P))*dx + delta1*dot(grad(q), grad(P))*dx + 0.001*P*q*dx
L = alpha1*dot(grad(q), grad(P0))*dx - q*div(uc)*dx - mult(delta1, dot( ugradu(uc, uc), grad(q))) *dx
#L = alpha1*dot(grad(q), grad(P0))*dx - q*div(uc)*dx
|