Numerical conductivity coefficient matrix

variational formulation
1

Hello,

I wish to create a numerical conductivity coefficient matrix, C, which would be part of the variational statement

a = dot(C*grad(u),grad(v))*dx

for TrialFunction u and TestFunction v.

I have the elements of the coefficient matrix at each node in the mesh stored in files and I intend to read the data into a numpy array, for example, for the numpy array c00:

c00 = read_data(“c00.dat”)
…

(read_data() is just some routine to get the numbers from the file c00.dat into the numpy array c00.)

I also know that I can create the matrix C using as_matrix:

C = as_matrix(((C00,C01,C02),(C01,C11,C12),(C02,C12,C22)))

(C is symmetric).

My question is how do I convert the numpy arrays c00, etc, into the quantities C00, etc, (and what type are the quantities C00)?

Perhaps I should not be using as_matrix in this instance? If so, what is the alternative? Can I simply write

C = ((c00,c01,c02),(c01,c11,c12),(c02,c12,c22))?

Also, should I have the conductivity values specified at each node or for each element in the mesh?

Thanks very much,

Peter.

1 Like
2

The following shows how to build components of the matrix as Function objects from the data. For P_1 functions the data are here given in mesh vertices while for P_0 they are given at cells.

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(4, 4)

# Fake data
x, y = mesh.coordinates().T
# Store and load for completness
np.savetxt('data.txt', np.c_[x, y])
c0, c1 = np.loadtxt('data.txt').T

# Want to have ci represent coefficient of P1 function
V = FunctionSpace(mesh, 'CG', 1)
# We have data ordered by vertex id while setting dof requires dof ordering
d2v = dof_to_vertex_map(V)

C0 = Function(V); C0.vector().set_local(c0[d2v])
C1 = Function(V); C1.vector().set_local(c1[d2v])

# Now C0 is like x, and C1 is like y
x, y = SpatialCoordinate(mesh)
assert assemble(inner(x-C0, x-C0)*dx) < 1E-13 and assemble(inner(y-C1, y-C1)*dx) < 1E-13

# Finally the matrix
C = as_matrix(((C0, 0), (0, C1)))

# P0 version
V = FunctionSpace(mesh, 'DG', 0)
# Data
x, y = (interpolate(Expression(e, degree=1), V).vector().get_local() for e in ('x[0]', 'x[1]'))
# Store and load for completness
np.savetxt('data.txt', np.c_[x, y])
c0, c1 = np.loadtxt('data.txt').T

# Dof to cell is here an identity, so d2c is c2d
c2d = V.dofmap().entity_dofs(mesh, mesh.topology().dim())
assert np.linalg.norm(c2d - np.arange(V.dim()), np.inf) < 1E-15
# Want to have ci represent coefficient of P0 function
C0 = Function(V); C0.vector().set_local(c0[c2d])
C1 = Function(V); C1.vector().set_local(c1[c2d])

# Now C0 is like x, and C1 is like y
x, y = SpatialCoordinate(mesh)
# With 0 order quadrature these should be same
dx0 = dx(metadata={'quadrature_degree': 0})
assert abs(assemble(inner(x-C0, x-C0)*dx0)) < 1E-13 and abs(assemble(inner(y-C1, y-C1)*dx0)) < 1E-13
3 Likes
3

Dear MiroK,

That is exactly what I wanted and it works very well.

Thank you very much.

Peter.