### How to assemble diagonal matrix from a Function-vector

269

views

0

Hello, everyone,

I want to assemble a diagonal matrix from a vector, but don't know how to code with python in FeniCS.

The question came from the enthalpy method dealing with solidification simulation, the assemmble equation is M{H}+K{T}=F, Here I declare three Fucntions : T, H, and dtdh; when solving with Newton Rapson method, the following equation comes, (M+K*[DTDH])dH = - R(H, T), In this equation the [DTDH] matrix is a diagnal matrix which is pointwisely calculated by T and H Function, My question is how to assemble DTDH matrix from dtdh function in FEniCS?

Before, I code the dtdh fucntion in the form to get the coeffcitiont matrix similar with K*[DTDH], however, dtdh here is declared as a Function, so this coeff matrix is not equal with K*[DTDH], so I want to change to the correct one.

Anyone konw how?

I want to assemble a diagonal matrix from a vector, but don't know how to code with python in FeniCS.

The question came from the enthalpy method dealing with solidification simulation, the assemmble equation is M{H}+K{T}=F, Here I declare three Fucntions : T, H, and dtdh; when solving with Newton Rapson method, the following equation comes, (M+K*[DTDH])dH = - R(H, T), In this equation the [DTDH] matrix is a diagnal matrix which is pointwisely calculated by T and H Function, My question is how to assemble DTDH matrix from dtdh function in FEniCS?

Before, I code the dtdh fucntion in the form to get the coeffcitiont matrix similar with K*[DTDH], however, dtdh here is declared as a Function, so this coeff matrix is not equal with K*[DTDH], so I want to change to the correct one.

Anyone konw how?

Community: FEniCS Project

### 1 Answer

3

Matrices have the method A.set_diagonal(vector). For the code you provided in your comment, the following should work. I didn't understand how you generate your dTdH function, but if you have one, you get a diagonal matrix like this:

```
from dolfin import *
mesh = UnitSquareMesh(5, 5)
V = FunctionSpace(mesh, 'P', 1)
du = TrialFunction(V)
v = TestFunction(V)
# dTdH = Function(V)
# set some dTdH values ...
dTdH = interpolate(Expression('x[0]*x[1]', degree=2), V)
M = assemble(du*v*dx)
M_dtdh = M.copy()
M_dtdh.zero()
M_dtdh.set_diagonal(dTdH.vector())
# check
x = dTdH.vector().copy()
x.zero()
M_dtdh.get_diagonal(x)
print(x.array())
```

Please login to add an answer/comment or follow this question.

"""

Simple version of the enthalpy method in solidification

H'= Laplace(T) + f in the unit square

u = u_D on the boundary

u = u_0 at t = 0

f = 1.0

Implicit integration is used here.

"""

from dolfin import *

# Create mesh and define function space

mesh = UnitSquareMesh(5, 5)

V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition

def boundary(x, on_boundary):

return on_boundary

bc = DirichletBC(V, 0.0, boundary)

# Define constant

alpha = Constant(1.0)

dt = Constant(0.1)

# Define initial value

H_n = interpolate(0.0, V)

# Define variational problem

du = TrialFunction(V)

v = TestFunction(V)

f = 1.0

# Define Functions to save current Temperature, Enthalpy and dTdH field

T = Function(V)

H = Function(V)

dTdH = Function(V)

# Update dTdH with the T and H nodal values

def UpdatedTdH(T, H, dTdH):

t_array = T.vector().array()

h_array = H.vector().array()

dtdh_array = dTdH.vector().array()

dtdh_array[:] = np.where(())

dTdH.vector()[:] = dtdh_array[:]

# Define a and l

# Method 1 The old way, just use dTdH function in the form, however not accurate

a = 1.0/dt*du*v*dx + dTdh*dot(grad(du), grad(v))*dx

l = 1.0/dt*(H-H_n)*v*dx + dot(grad(T), grad(v))*dx-f*v*dx

# Method 2 The new way, assembly matrix M, K and dTdH to form the lhs matrix

M = assemble(du*v*dx)

K = assemble(dot(grad(du), grad(v))*dx)

# Then assemble the dTdH diagonal matrix with the dTdH Function nodal value

M_dtdh=??????

# Finally the lhs matrix, Then solve .

A = M/dt+K*M_dTdH