# Best way to specify coefficient fields

Asked by Chris Richardson on 2010-02-26

What is the best way to specify coefficient fields, e.g. conductivity, permeability, elastic constants ?
I have seen the example "tensor-weighted-poisson" which seems very convoluted.

I am interested in solving Darcy flow with variable permeability k (usually a scalar, but may be a tensor with diagonal entries)...
Also, elastic problems with anisotropic material whose properties vary through space.
There are 9 independent coefficients for elasticity, (out of 81 entries in Cijkl, which relates stress s_ij and strain e_kl).

## Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Chris Richardson
Solved:
2010-02-27
Last query:
2010-02-27
2010-02-26
 Anders Logg (logg) said on 2010-02-26: #1

On Fri, Feb 26, 2010 at 11:45:10AM -0000, Chris Richardson wrote:
> New question #102501 on DOLFIN:
>
> What is the best way to specify coefficient fields, e.g. conductivity, permeability, elastic constants ?
> I have seen the example "tensor-weighted-poisson" which seems very convoluted.
>
> I am interested in solving Darcy flow with variable permeability k (usually a scalar, but may be a tensor with diagonal entries)...
> Also, elastic problems with anisotropic material whose properties vary through space.
> There are 9 independent coefficients for elasticity, (out of 81 entries in Cijkl, which relates stress s_ij and strain e_kl).

It depends on where your data for the field comes from.

The reason for the tensor-weighted-poisson demo being a bit complex is
that in addition to demonstrating tensor-valued coefficients, it also
demonstrates

1. Use of C++ code from Python (for efficiency)
2. Reading cell-based field data from file

--
Anders

 Chris Richardson (chris-bpi) said on 2010-02-26: #2

I would like to have the coefficient values for each point in space stored in memory, somehow.
I suppose the values might be originally read from a file, but then might be altered during the course of the program...

Sorry if it seems a basic question, I am only just starting out with Dolfin...

 Harish Narayanan (hnarayanan) said on 2010-02-26: #3

Will this work for you? e.g, for inverse of the permeability

k = "1.0/(exp(-(((x[1] - 0.5 - 0.1*sin(10*x[0]))/0.1)*((x[1] - 0.5 - 0.1*sin(10*x[0]))/0.1))) + 1.0)"
kinv11 = Expression(k)
kinv12 = Constant(0.0)
kinv21 = Constant(0.0)
kinv22 = Expression(k)
Kinv = as_matrix(((kinv11, kinv12), (kinv21, kinv22)))

Of course, this needn't be isotropic.

And later on, when defining your weak form, you'll have a term like the following:

a = inner(v, Kinv*u)*dx

 Chris Richardson (chris-bpi) said on 2010-02-26: #4

Not really. I want k to be a variable, not an expression.

 Garth Wells (garth-wells) said on 2010-02-26: #5

On Fri, 2010-02-26 at 13:06 +0000, Chris Richardson wrote:
> Question #102501 on DOLFIN changed:
>
> Chris Richardson posted a new comment:
> Not really. I want k to be a variable, not an expression.
>

Expressions can vary- you can make them a function of space, time, etc.
If you want something to be stored in a Vector, you can project it onto
a finite element space or interpolate it in a FE space. For example, if
you project onto a P0 (zeroth order) element, you'll have something
which varies cell-wise.

Garth

 Chris Richardson (chris-bpi) said on 2010-02-27: #6

OK, so I can use something like this for a scalar field defined on the nodes:

V = FunctionSpace(mesh, "CG", 1)
# create an initial condition for C
initC = Expression("0.5*tanh((x[1]-1)*20)+1.5")
C = project(initC,V)

and presumably, the same will work with Vector and Tensor Functionspaces...

Now, I might have some other questions, but I think this is solved for me...