On Wed, 2006-02-22 at 10:45 +0000, Alexander Jarosch wrote:
Hi,
I have a rather simple question I guess. I tried to figure out how to
give the stokes solver a viscosity for each vertex and started of with
the stokes solver with a constant viscosity nu:
modified from Anders file:
---------------
P1 = FiniteElement("Lagrange", "triangle", 1)
P2 = FiniteElement("Vector Lagrange", "triangle", 2)
TH = P2 + P1
(v, q) = BasisFunctions(TH)
(u, p) = BasisFunctions(TH)
nu = Constant() #viscosity
f = Function(P2)
a = (dot(nu*grad(u), grad(v)) - p*div(v) + div(u)*q)*dx
L = dot(f, v)*dx
---------------
now that works fine and now i have the viscosity values for each vertex
in a vector as values and in a function like:
Function viscfunc(viscosity, mesh, element);
One other thing, this is implicitly a DiscreteFunction. Trying making a
constant as I suggested, and as a next step try a user-defined Function
to specify nu as a function of position, in the same way the source term
is defined in src/demo/solvers/stokes.
Garth
Now i tried to get a form file which takes these values like:
---------------
P1 = FiniteElement("Lagrange", "triangle", 1)
P2 = FiniteElement("Vector Lagrange", "triangle", 2)
element1 = FiniteElement("Discontinuous vector Lagrange", "triangle", 0, 1)
TH = P2 + P1
(v, q) = BasisFunctions(TH)
(u, p) = BasisFunctions(TH)
nu = Function(element1) #viscosity
f = Function(P2)
a = (dot(dot(nu, grad(u)), grad(v)) - p*div(v) + div(u)*q)*dx
L = dot(f, v)*dx
---------------
well and it does not work, as well as the approach like:
---------------
P1 = FiniteElement("Lagrange", "triangle", 1)
P2 = FiniteElement("Vector Lagrange", "triangle", 2)
constant_scalar = FiniteElement("Discontinuous Lagrange", "triangle", 1)
TH = P2 + P1
(v, q) = BasisFunctions(TH)
(u, p) = BasisFunctions(TH)
nu = Function(constant_scalar) #viscosity
#nu = Constant()
f = Function(P2)
a = (dot(nu*grad(u), grad(v)) - p*div(v) + div(u)*q)*dx
L = dot(f, v)*dx
---------------
which actually takes my viscosity vector and starts to solve but does
not converge.
might be a simple mistake in there, I appreciate any help or hints.
Alex
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
