7
$\begingroup$

I am continuing my dive into computational fluid dynamics. I would like to build a simple test case for modelling the flow around a cylinder in FEniCS and then continue to model turbulence in order to observe the Karman vortex street.

I am a newby in both FEniCS and CFD. I have started to go through the FEniCS book and went through the first simple scripts in the tutorial.

Since this is a standard problem and I do not like to re-invent the wheel, my question is: is there some tutorial available on how to do this? Even a not-well commented code / FEniCS script would do as a reference.

I have googled this and did not get much out of it, which is why I ask here.

Improve this question
$\endgroup$

3 Answers 3

Reset to default
11
$\begingroup$

I suggest you start by looking at the FEniCS Navier-Stokes demo which is documented here:

http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/navier-stokes/python/documentation.html

For your specific test case, you might want to look at the NSBENCH set of test problem (which are described in Chapter 21 of the FEniCS Book). The code for that book chapter is available at

http://launchpad.net/nsbench/

The code for your test case (flow around a cylinder) can be found here:

http://bazaar.launchpad.net/~nsbench/nsbench/main/view/head:/problems/cylinder.py

Note that this code depends on the interface defined by the solvers in NSBENCH, but it might still be useful to look at for how to set the boundary conditions, and how to design an interface for your own solver.

Improve this answer
$\endgroup$
5
$\begingroup$

This is I think rotating cylinder and Karman vortex street behind it. You can also uncomment SUPG/PSPG stabilization. But I suspect that mesh size h computed as 2*circumradius is not good quantity and prefer cell diameter computed this way https://scicomp.stackexchange.com/a/7181/4254

from dolfin import *
import time

#set_log_level(PROGRESS)
parameters['form_compiler']['representation'] = 'quadrature'
parameters['form_compiler']['optimize'] = True
parameters['form_compiler']['cpp_optimize'] = True

center = Point(0.2, 0.2)
radius = 0.05

class Cylinder(SubDomain):
    def inside(self,x,on_boundary):
        r = sqrt((x[0] - center[0])**2 + (x[1] - center[1])**2)
        return(on_boundary and (r < 2*radius + sqrt(DOLFIN_EPS)))
    def snap(self, x):
        r = sqrt((x[0] - center[0])**2 + (x[1] - center[1])**2)
        if r <= radius:
            x[0] = center[0] + (radius / r)*(x[0] - center[0])
            x[1] = center[1] + (radius / r)*(x[1] - center[1])

class Walls(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and ((x[1] < DOLFIN_EPS) or (x[1] > (0.41 - DOLFIN_EPS))))

class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and (x[0] < DOLFIN_EPS))

class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and (x[0] > (2.2 - DOLFIN_EPS)))

# Parameters
nu = Constant(0.001)
dt = 0.1
idt = Constant(1./dt)
t_end = 10.0
theta=0.5   # Crank-Nicholson timestepping

# Mesh
mesh = Mesh("bench_L1.xml")
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(0)
Gamma = Cylinder()
Gamma.mark(boundary_parts, 1)

#refine the mesh and snap boundary
mesh = refine(mesh)
mesh.snap_boundary(Gamma, False)

# Define function spaces (Taylor-Hood)
V = VectorFunctionSpace(mesh, "CG", 2)
P = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V, P])

# No-slip boundary condition for velocity
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, Walls())

vc = Expression(("-cos(atan2(x[0]-0.2,x[1]-0.2))","sin(atan2(x[0]-0.2,x[1]-0.2))"))
bc_cylinder = DirichletBC(W.sub(0), vc, Cylinder())

# Inflow boundary condition for velocity and temperature
v_in = Expression(("1.5 * 4.0 * x[1] * (0.41 - x[1]) / ( 0.41 * 0.41 )","0.0"))
bc1 = DirichletBC(W.sub(0), v_in, Inflow())

# Collect boundary conditions
bcs = [bc_cylinder, bc0, bc1]

# Define unknown and test function(s)
(u, q) = TestFunctions(W)

n = FacetNormal(mesh)
I = Identity(V.cell().d)    # Identity tensor

# current time step
w = Function(W)
(v, p) = (as_vector((w[0], w[1])), w[2])
D = 0.5*(grad(v)+grad(v).T)
T = -p*I + 2.0*nu*D

# previous time step
w0 = Function(W)
(v0, p0) = (as_vector((w0[0], w0[1])), w0[2])
D0 = 0.5*(grad(v0)+grad(v0).T)
T0 = -p*I + 2.0*nu*D0

# Define variational forms without time derivative in previous time
F0_eq1 = (inner(T0, grad(u)))*dx + inner(grad(v0)*v0, u)*dx
F0_eq2 = 0*q*dx
F0 = F0_eq1 + F0_eq2

# variational form without time derivative in current time
F1_eq1 = (inner(T, grad(u)) + inner(grad(v)*v, u))*dx
F1_eq2 = q*div(v)*dx
F1 = F1_eq1 + F1_eq2

# combine variational forms with time derivative
#
#  dw/dt + F(t) = 0 is approximated as
#  (w-w0)/dt + (1-theta)*F(t0) + theta*F(t) = 0
#
F = idt*inner((v-v0),u)*dx + (1.0-theta)*F0 + theta*F1

# residual of strong Navier-Stokes
r = idt*(v-v0) + theta*grad(v)*v + (1.0-theta)*grad(v0)*v0 \
   - theta*div(T) - (1.0-theta)*div(T0)

# stabilization parameter
h = CellSize(mesh)
velocity = v0
vnorm = sqrt(dot(velocity, velocity))
tau = ( (2.0*idt)**2 + (2.0*vnorm/h)**2 + (4.0*nu/h**2)**2 )**(-0.5)

# add SUPG stabilization
#F += tau*inner(grad(u)*v, r)*dx

# add PSPG stabilization
#F += tau*inner(grad(q), r)*dx

# define Jacobian
J = derivative(F, w)

# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")

# create variational problem and solver
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['maximum_iterations'] = 20

# Time-stepping
t = dt
while t < t_end:

    print "t =", t

    # Compute
    begin("Solving ....")
    solver.solve()
    end()

    # Extract solutions:
    (v, p) = w.split()

    # Plot
    #plot(v)

    # Save to file
    ufile << v
    pfile << p

    # Move to next time step
    w0.assign(w)
    t += dt

and mesh file

<?xml version="1.0"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
  <mesh celltype="triangle" dim="2">
    <vertices size="147">
      <vertex index="0" x="0" y="0" />
      <vertex index="1" x="2.2" y="0" />
      <vertex index="2" x="2.2" y="0.41" />
      <vertex index="3" x="0" y="0.41" />
      <vertex index="4" x="0.15" y="0.2" />
      <vertex index="5" x="0.25" y="0.2" />
      <vertex index="6" x="0.2" y="0.15" />
      <vertex index="7" x="0.2" y="0.25" />
      <vertex index="8" x="0.104762" y="0" />
      <vertex index="9" x="0.209524" y="0" />
      <vertex index="10" x="0.314286" y="0" />
      <vertex index="11" x="0.419048" y="0" />
      <vertex index="12" x="0.52381" y="0" />
      <vertex index="13" x="0.628571" y="0" />
      <vertex index="14" x="0.733333" y="0" />
      <vertex index="15" x="0.838095" y="0" />
      <vertex index="16" x="0.942857" y="0" />
      <vertex index="17" x="1.04762" y="0" />
      <vertex index="18" x="1.15238" y="0" />
      <vertex index="19" x="1.25714" y="0" />
      <vertex index="20" x="1.3619" y="0" />
      <vertex index="21" x="1.46667" y="0" />
      <vertex index="22" x="1.57143" y="0" />
      <vertex index="23" x="1.67619" y="0" />
      <vertex index="24" x="1.78095" y="0" />
      <vertex index="25" x="1.88571" y="0" />
      <vertex index="26" x="1.99048" y="0" />
      <vertex index="27" x="2.09524" y="0" />
      <vertex index="28" x="2.2" y="0.1025" />
      <vertex index="29" x="2.2" y="0.205" />
      <vertex index="30" x="2.2" y="0.3075" />
      <vertex index="31" x="2.1" y="0.41" />
      <vertex index="32" x="2" y="0.41" />
      <vertex index="33" x="1.9" y="0.41" />
      <vertex index="34" x="1.8" y="0.41" />
      <vertex index="35" x="1.7" y="0.41" />
      <vertex index="36" x="1.6" y="0.41" />
      <vertex index="37" x="1.5" y="0.41" />
      <vertex index="38" x="1.4" y="0.41" />
      <vertex index="39" x="1.3" y="0.41" />
      <vertex index="40" x="1.2" y="0.41" />
      <vertex index="41" x="1.1" y="0.41" />
      <vertex index="42" x="1" y="0.41" />
      <vertex index="43" x="0.9" y="0.41" />
      <vertex index="44" x="0.8" y="0.41" />
      <vertex index="45" x="0.7" y="0.41" />
      <vertex index="46" x="0.6" y="0.41" />
      <vertex index="47" x="0.5" y="0.41" />
      <vertex index="48" x="0.4" y="0.41" />
      <vertex index="49" x="0.3" y="0.41" />
      <vertex index="50" x="0.2" y="0.41" />
      <vertex index="51" x="0.1" y="0.41" />
      <vertex index="52" x="0" y="0.3075" />
      <vertex index="53" x="0" y="0.205" />
      <vertex index="54" x="0" y="0.1025" />
      <vertex index="55" x="0.156699" y="0.225" />
      <vertex index="56" x="0.175" y="0.243301" />
      <vertex index="57" x="0.225" y="0.243301" />
      <vertex index="58" x="0.243301" y="0.225" />
      <vertex index="59" x="0.243301" y="0.175" />
      <vertex index="60" x="0.225" y="0.156699" />
      <vertex index="61" x="0.175" y="0.156699" />
      <vertex index="62" x="0.156699" y="0.175" />
      <vertex index="63" x="0.352033" y="0.262216" />
      <vertex index="64" x="1.04579" y="0.252253" />
      <vertex index="65" x="1.24065" y="0.243419" />
      <vertex index="66" x="0.985775" y="0.171509" />
      <vertex index="67" x="1.30324" y="0.187395" />
      <vertex index="68" x="0.890988" y="0.190768" />
      <vertex index="69" x="1.39085" y="0.202767" />
      <vertex index="70" x="0.801464" y="0.204685" />
      <vertex index="71" x="1.49904" y="0.204666" />
      <vertex index="72" x="0.704267" y="0.202289" />
      <vertex index="73" x="1.60238" y="0.204951" />
      <vertex index="74" x="0.487789" y="0.19191" />
      <vertex index="75" x="1.72935" y="0.219637" />
      <vertex index="76" x="1.83335" y="0.187427" />
      <vertex index="77" x="2.03359" y="0.193642" />
      <vertex index="78" x="0.314663" y="0.307348" />
      <vertex index="79" x="0.249039" y="0.333583" />
      <vertex index="80" x="2.12147" y="0.253405" />
      <vertex index="81" x="2.06373" y="0.0836837" />
      <vertex index="82" x="0.0734726" y="0.267055" />
      <vertex index="83" x="0.143989" y="0.0677174" />
      <vertex index="84" x="0.0613613" y="0.33781" />
      <vertex index="85" x="0.299628" y="0.236829" />
      <vertex index="86" x="0.527114" y="0.295386" />
      <vertex index="87" x="0.262384" y="0.26549" />
      <vertex index="88" x="0.275777" y="0.18341" />
      <vertex index="89" x="1.74578" y="0.330928" />
      <vertex index="90" x="1.84001" y="0.302448" />
      <vertex index="91" x="1.64998" y="0.313838" />
      <vertex index="92" x="0.648421" y="0.302367" />
      <vertex index="93" x="1.54544" y="0.308508" />
      <vertex index="94" x="0.755686" y="0.306156" />
      <vertex index="95" x="1.42668" y="0.314809" />
      <vertex index="96" x="0.870524" y="0.298596" />
      <vertex index="97" x="1.22553" y="0.112785" />
      <vertex index="98" x="1.05603" y="0.147372" />
      <vertex index="99" x="0.969346" y="0.329867" />
      <vertex index="100" x="1.32527" y="0.276284" />
      <vertex index="101" x="2.02914" y="0.29511" />
      <vertex index="102" x="1.88168" y="0.119486" />
      <vertex index="103" x="1.79545" y="0.104842" />
      <vertex index="104" x="1.35688" y="0.0892539" />
      <vertex index="105" x="1.69041" y="0.105324" />
      <vertex index="106" x="0.914199" y="0.0955566" />
      <vertex index="107" x="1.58234" y="0.102923" />
      <vertex index="108" x="0.61843" y="0.10815" />
      <vertex index="109" x="1.47999" y="0.101161" />
      <vertex index="110" x="0.725217" y="0.104273" />
      <vertex index="111" x="0.821727" y="0.102306" />
      <vertex index="112" x="1.96725" y="0.105778" />
      <vertex index="113" x="0.418624" y="0.241547" />
      <vertex index="114" x="0.135226" y="0.239087" />
      <vertex index="115" x="1.14921" y="0.314841" />
      <vertex index="116" x="0.404508" y="0.333259" />
      <vertex index="117" x="0.590155" y="0.206269" />
      <vertex index="118" x="2.12554" y="0.150731" />
      <vertex index="119" x="1.93166" y="0.22436" />
      <vertex index="120" x="0.104689" y="0.190046" />
      <vertex index="121" x="0.267856" y="0.0744753" />
      <vertex index="122" x="0.269368" y="0.221581" />
      <vertex index="123" x="0.434446" y="0.114927" />
      <vertex index="124" x="1.0089" y="0.0793381" />
      <vertex index="125" x="0.953023" y="0.246698" />
      <vertex index="126" x="1.25853" y="0.33355" />
      <vertex index="127" x="1.05268" y="0.340828" />
      <vertex index="128" x="2.1369" y="0.067383" />
      <vertex index="129" x="0.522398" y="0.112537" />
      <vertex index="130" x="1.3421" y="0.348929" />
      <vertex index="131" x="0.215186" y="0.289473" />
      <vertex index="132" x="0.203197" y="0.335348" />
      <vertex index="133" x="1.13984" y="0.190728" />
      <vertex index="134" x="1.93813" y="0.330504" />
      <vertex index="135" x="1.10505" y="0.0883705" />
      <vertex index="136" x="0.341988" y="0.0824728" />
      <vertex index="137" x="0.253113" y="0.138035" />
      <vertex index="138" x="0.291932" y="0.131907" />
      <vertex index="139" x="0.203713" y="0.102439" />
      <vertex index="140" x="0.150375" y="0.134668" />
      <vertex index="141" x="0.151494" y="0.309119" />
      <vertex index="142" x="0.354384" y="0.17896" />
      <vertex index="143" x="0.0685367" y="0.100776" />
      <vertex index="144" x="0.424128" y="0.180088" />
      <vertex index="145" x="0.474925" y="0.0568661" />
      <vertex index="146" x="2.10843" y="0.347669" />
    </vertices>
    <cells size="231">
      <triangle index="0" v0="29" v1="30" v2="80" />
      <triangle index="1" v0="26" v1="27" v2="81" />
      <triangle index="2" v0="8" v1="9" v2="83" />
      <triangle index="3" v0="3" v1="51" v2="84" />
      <triangle index="4" v0="3" v1="52" v2="84" />
      <triangle index="5" v0="49" v1="78" v2="79" />
      <triangle index="6" v0="52" v1="82" v2="84" />
      <triangle index="7" v0="63" v1="78" v2="85" />
      <triangle index="8" v0="78" v1="79" v2="87" />
      <triangle index="9" v0="77" v1="80" v2="101" />
      <triangle index="10" v0="26" v1="81" v2="112" />
      <triangle index="11" v0="77" v1="81" v2="112" />
      <triangle index="12" v0="74" v1="86" v2="113" />
      <triangle index="13" v0="28" v1="29" v2="118" />
      <triangle index="14" v0="76" v1="102" v2="119" />
      <triangle index="15" v0="76" v1="90" v2="119" />
      <triangle index="16" v0="72" v1="108" v2="117" />
      <triangle index="17" v0="63" v1="78" v2="116" />
      <triangle index="18" v0="55" v1="56" v2="114" />
      <triangle index="19" v0="5" v1="59" v2="88" />
      <triangle index="20" v0="78" v1="85" v2="87" />
      <triangle index="21" v0="64" v1="66" v2="98" />
      <triangle index="22" v0="65" v1="67" v2="97" />
      <triangle index="23" v0="67" v1="97" v2="104" />
      <triangle index="24" v0="15" v1="16" v2="106" />
      <triangle index="25" v0="66" v1="68" v2="106" />
      <triangle index="26" v0="68" v1="106" v2="111" />
      <triangle index="27" v0="15" v1="106" v2="111" />
      <triangle index="28" v0="20" v1="21" v2="104" />
      <triangle index="29" v0="67" v1="69" v2="104" />
      <triangle index="30" v0="69" v1="104" v2="109" />
      <triangle index="31" v0="21" v1="104" v2="109" />
      <triangle index="32" v0="14" v1="15" v2="111" />
      <triangle index="33" v0="68" v1="70" v2="111" />
      <triangle index="34" v0="70" v1="110" v2="111" />
      <triangle index="35" v0="14" v1="110" v2="111" />
      <triangle index="36" v0="21" v1="22" v2="109" />
      <triangle index="37" v0="69" v1="71" v2="109" />
      <triangle index="38" v0="71" v1="107" v2="109" />
      <triangle index="39" v0="22" v1="107" v2="109" />
      <triangle index="40" v0="77" v1="101" v2="119" />
      <triangle index="41" v0="33" v1="34" v2="90" />
      <triangle index="42" v0="75" v1="76" v2="90" />
      <triangle index="43" v0="34" v1="35" v2="89" />
      <triangle index="44" v0="34" v1="89" v2="90" />
      <triangle index="45" v0="75" v1="89" v2="90" />
      <triangle index="46" v0="35" v1="36" v2="91" />
      <triangle index="47" v0="35" v1="89" v2="91" />
      <triangle index="48" v0="73" v1="75" v2="91" />
      <triangle index="49" v0="75" v1="89" v2="91" />
      <triangle index="50" v0="45" v1="46" v2="92" />
      <triangle index="51" v0="72" v1="92" v2="117" />
      <triangle index="52" v0="36" v1="37" v2="93" />
      <triangle index="53" v0="36" v1="91" v2="93" />
      <triangle index="54" v0="71" v1="73" v2="93" />
      <triangle index="55" v0="73" v1="91" v2="93" />
      <triangle index="56" v0="44" v1="45" v2="94" />
      <triangle index="57" v0="45" v1="92" v2="94" />
      <triangle index="58" v0="70" v1="72" v2="94" />
      <triangle index="59" v0="72" v1="92" v2="94" />
      <triangle index="60" v0="37" v1="38" v2="95" />
      <triangle index="61" v0="37" v1="93" v2="95" />
      <triangle index="62" v0="69" v1="71" v2="95" />
      <triangle index="63" v0="71" v1="93" v2="95" />
      <triangle index="64" v0="43" v1="44" v2="96" />
      <triangle index="65" v0="44" v1="94" v2="96" />
      <triangle index="66" v0="68" v1="70" v2="96" />
      <triangle index="67" v0="70" v1="94" v2="96" />
      <triangle index="68" v0="42" v1="43" v2="99" />
      <triangle index="69" v0="43" v1="96" v2="99" />
      <triangle index="70" v0="13" v1="14" v2="110" />
      <triangle index="71" v0="70" v1="72" v2="110" />
      <triangle index="72" v0="72" v1="108" v2="110" />
      <triangle index="73" v0="13" v1="108" v2="110" />
      <triangle index="74" v0="22" v1="23" v2="107" />
      <triangle index="75" v0="71" v1="73" v2="107" />
      <triangle index="76" v0="73" v1="105" v2="107" />
      <triangle index="77" v0="23" v1="105" v2="107" />
      <triangle index="78" v0="23" v1="24" v2="105" />
      <triangle index="79" v0="73" v1="75" v2="105" />
      <triangle index="80" v0="75" v1="103" v2="105" />
      <triangle index="81" v0="24" v1="103" v2="105" />
      <triangle index="82" v0="24" v1="25" v2="103" />
      <triangle index="83" v0="75" v1="76" v2="103" />
      <triangle index="84" v0="76" v1="102" v2="103" />
      <triangle index="85" v0="25" v1="102" v2="103" />
      <triangle index="86" v0="40" v1="41" v2="115" />
      <triangle index="87" v0="65" v1="67" v2="100" />
      <triangle index="88" v0="67" v1="69" v2="100" />
      <triangle index="89" v0="69" v1="95" v2="100" />
      <triangle index="90" v0="25" v1="26" v2="112" />
      <triangle index="91" v0="25" v1="102" v2="112" />
      <triangle index="92" v0="77" v1="112" v2="119" />
      <triangle index="93" v0="102" v1="112" v2="119" />
      <triangle index="94" v0="48" v1="49" v2="116" />
      <triangle index="95" v0="49" v1="78" v2="116" />
      <triangle index="96" v0="46" v1="47" v2="86" />
      <triangle index="97" v0="46" v1="86" v2="92" />
      <triangle index="98" v0="74" v1="86" v2="117" />
      <triangle index="99" v0="86" v1="92" v2="117" />
      <triangle index="100" v0="47" v1="48" v2="116" />
      <triangle index="101" v0="47" v1="86" v2="116" />
      <triangle index="102" v0="63" v1="113" v2="116" />
      <triangle index="103" v0="86" v1="113" v2="116" />
      <triangle index="104" v0="82" v1="114" v2="120" />
      <triangle index="105" v0="49" v1="50" v2="79" />
      <triangle index="106" v0="9" v1="10" v2="121" />
      <triangle index="107" v0="58" v1="87" v2="122" />
      <triangle index="108" v0="85" v1="87" v2="122" />
      <triangle index="109" v0="66" v1="98" v2="124" />
      <triangle index="110" v0="66" v1="106" v2="124" />
      <triangle index="111" v0="66" v1="68" v2="125" />
      <triangle index="112" v0="96" v1="99" v2="125" />
      <triangle index="113" v0="68" v1="96" v2="125" />
      <triangle index="114" v0="39" v1="40" v2="126" />
      <triangle index="115" v0="41" v1="42" v2="127" />
      <triangle index="116" v0="42" v1="99" v2="127" />
      <triangle index="117" v0="1" v1="27" v2="128" />
      <triangle index="118" v0="1" v1="28" v2="128" />
      <triangle index="119" v0="27" v1="81" v2="128" />
      <triangle index="120" v0="108" v1="117" v2="129" />
      <triangle index="121" v0="74" v1="117" v2="129" />
      <triangle index="122" v0="38" v1="39" v2="130" />
      <triangle index="123" v0="95" v1="100" v2="130" />
      <triangle index="124" v0="38" v1="95" v2="130" />
      <triangle index="125" v0="64" v1="66" v2="125" />
      <triangle index="126" v0="64" v1="99" v2="125" />
      <triangle index="127" v0="39" v1="126" v2="130" />
      <triangle index="128" v0="100" v1="126" v2="130" />
      <triangle index="129" v0="28" v1="118" v2="128" />
      <triangle index="130" v0="81" v1="118" v2="128" />
      <triangle index="131" v0="5" v1="88" v2="122" />
      <triangle index="132" v0="5" v1="58" v2="122" />
      <triangle index="133" v0="85" v1="88" v2="122" />
      <triangle index="134" v0="29" v1="80" v2="118" />
      <triangle index="135" v0="77" v1="81" v2="118" />
      <triangle index="136" v0="77" v1="80" v2="118" />
      <triangle index="137" v0="16" v1="17" v2="124" />
      <triangle index="138" v0="16" v1="106" v2="124" />
      <triangle index="139" v0="64" v1="99" v2="127" />
      <triangle index="140" v0="41" v1="115" v2="127" />
      <triangle index="141" v0="64" v1="115" v2="127" />
      <triangle index="142" v0="12" v1="13" v2="108" />
      <triangle index="143" v0="12" v1="108" v2="129" />
      <triangle index="144" v0="65" v1="100" v2="126" />
      <triangle index="145" v0="65" v1="115" v2="126" />
      <triangle index="146" v0="40" v1="115" v2="126" />
      <triangle index="147" v0="79" v1="87" v2="131" />
      <triangle index="148" v0="50" v1="79" v2="132" />
      <triangle index="149" v0="65" v1="115" v2="133" />
      <triangle index="150" v0="64" v1="115" v2="133" />
      <triangle index="151" v0="79" v1="131" v2="132" />
      <triangle index="152" v0="57" v1="87" v2="131" />
      <triangle index="153" v0="57" v1="58" v2="87" />
      <triangle index="154" v0="74" v1="123" v2="129" />
      <triangle index="155" v0="32" v1="33" v2="134" />
      <triangle index="156" v0="32" v1="101" v2="134" />
      <triangle index="157" v0="33" v1="90" v2="134" />
      <triangle index="158" v0="90" v1="119" v2="134" />
      <triangle index="159" v0="101" v1="119" v2="134" />
      <triangle index="160" v0="65" v1="97" v2="133" />
      <triangle index="161" v0="64" v1="98" v2="133" />
      <triangle index="162" v0="17" v1="18" v2="135" />
      <triangle index="163" v0="17" v1="124" v2="135" />
      <triangle index="164" v0="98" v1="124" v2="135" />
      <triangle index="165" v0="97" v1="133" v2="135" />
      <triangle index="166" v0="98" v1="133" v2="135" />
      <triangle index="167" v0="10" v1="121" v2="136" />
      <triangle index="168" v0="10" v1="11" v2="136" />
      <triangle index="169" v0="11" v1="123" v2="136" />
      <triangle index="170" v0="121" v1="136" v2="138" />
      <triangle index="171" v0="121" v1="137" v2="138" />
      <triangle index="172" v0="60" v1="137" v2="139" />
      <triangle index="173" v0="19" v1="20" v2="104" />
      <triangle index="174" v0="19" v1="97" v2="104" />
      <triangle index="175" v0="121" v1="137" v2="139" />
      <triangle index="176" v0="18" v1="19" v2="97" />
      <triangle index="177" v0="18" v1="97" v2="135" />
      <triangle index="178" v0="9" v1="121" v2="139" />
      <triangle index="179" v0="9" v1="83" v2="139" />
      <triangle index="180" v0="52" v1="53" v2="82" />
      <triangle index="181" v0="59" v1="88" v2="137" />
      <triangle index="182" v0="59" v1="60" v2="137" />
      <triangle index="183" v0="88" v1="137" v2="138" />
      <triangle index="184" v0="4" v1="62" v2="120" />
      <triangle index="185" v0="62" v1="120" v2="140" />
      <triangle index="186" v0="53" v1="82" v2="120" />
      <triangle index="187" v0="82" v1="114" v2="141" />
      <triangle index="188" v0="85" v1="88" v2="142" />
      <triangle index="189" v0="88" v1="138" v2="142" />
      <triangle index="190" v0="0" v1="8" v2="143" />
      <triangle index="191" v0="0" v1="54" v2="143" />
      <triangle index="192" v0="74" v1="113" v2="144" />
      <triangle index="193" v0="74" v1="123" v2="144" />
      <triangle index="194" v0="82" v1="84" v2="141" />
      <triangle index="195" v0="63" v1="85" v2="142" />
      <triangle index="196" v0="8" v1="83" v2="143" />
      <triangle index="197" v0="53" v1="54" v2="143" />
      <triangle index="198" v0="83" v1="140" v2="143" />
      <triangle index="199" v0="120" v1="140" v2="143" />
      <triangle index="200" v0="53" v1="120" v2="143" />
      <triangle index="201" v0="61" v1="62" v2="140" />
      <triangle index="202" v0="136" v1="138" v2="142" />
      <triangle index="203" v0="63" v1="113" v2="142" />
      <triangle index="204" v0="113" v1="142" v2="144" />
      <triangle index="205" v0="123" v1="136" v2="142" />
      <triangle index="206" v0="123" v1="142" v2="144" />
      <triangle index="207" v0="50" v1="51" v2="141" />
      <triangle index="208" v0="51" v1="84" v2="141" />
      <triangle index="209" v0="50" v1="132" v2="141" />
      <triangle index="210" v0="4" v1="55" v2="114" />
      <triangle index="211" v0="4" v1="114" v2="120" />
      <triangle index="212" v0="11" v1="12" v2="145" />
      <triangle index="213" v0="123" v1="129" v2="145" />
      <triangle index="214" v0="11" v1="123" v2="145" />
      <triangle index="215" v0="12" v1="129" v2="145" />
      <triangle index="216" v0="6" v1="60" v2="139" />
      <triangle index="217" v0="131" v1="132" v2="141" />
      <triangle index="218" v0="6" v1="61" v2="140" />
      <triangle index="219" v0="83" v1="139" v2="140" />
      <triangle index="220" v0="6" v1="139" v2="140" />
      <triangle index="221" v0="2" v1="30" v2="146" />
      <triangle index="222" v0="2" v1="31" v2="146" />
      <triangle index="223" v0="31" v1="32" v2="146" />
      <triangle index="224" v0="80" v1="101" v2="146" />
      <triangle index="225" v0="30" v1="80" v2="146" />
      <triangle index="226" v0="32" v1="101" v2="146" />
      <triangle index="227" v0="7" v1="57" v2="131" />
      <triangle index="228" v0="56" v1="114" v2="141" />
      <triangle index="229" v0="7" v1="56" v2="141" />
      <triangle index="230" v0="7" v1="131" v2="141" />
    </cells>
  </mesh>
</dolfin>
Improve this answer
$\endgroup$
3
  • $\begingroup$ You can also add interior penalty stabilization as described in scicomp.stackexchange.com/a/7181/4254 . It seems more robust with this code than SUPG/PSPG as implemented here. $\endgroup$
    – Jan Blechta
    May 20, 2013 at 19:26
  • $\begingroup$ when I try running your code I get: RuntimeError: Cannot create a MeshFunction of type 'size_t'. Since I'm a complete FEniCS newby, I don't know wha to make of this. $\endgroup$
    – seb
    May 22, 2013 at 6:28
  • $\begingroup$ @sebastian: In newer versions of FEniCS size_t type is used in place of uint type. Replace size_t for uint. $\endgroup$
    – Jan Blechta
    May 22, 2013 at 10:31
-2
$\begingroup$

There is a problem with a SUPG / PSPG implementation. The test does not converge. How do I solve this problem?

Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.