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

269
views
0
9 months ago by
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?

Community: FEniCS Project
written 9 months ago by Miguel
Sorry, The ref code is added here, can't run now, but may make my question clear

"""
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

# Method 2 The new way, assembly matrix M, K and dTdH to form the lhs matrix
M = assemble(du*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
written 9 months ago by neochen

3
9 months ago by
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())​