This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Efficient assembling of a matrix out of an existing matrix

+2 votes

I am trying to implement a FEM-FCT algorithm for solving a transport equation for a special purpose. Source of the algorithm: http://www.mathematik.uni-dortmund.de/~kuzmin/linfct.pdf

My specific problem is, that I have a known Matrix $K_{ij}$ and I have to assemble a new matrix $D_{ij}$ using components of $K_{ij}$ the following way:

$$D_{ij}=max(-K_{ij},0,-K_{ji})$$

  1. Is there a efficient/easy way to do this?

  2. I have found a dirty way to achieve it, using the set command. But I guess it is slow and it also seems that there is a bug:
    I write blocks of components into the new matrix. It is working for a 3x3 block, while it produces an error for others, like 4x4. This seems to me like it is not working as intended? (See code below.)

from fenics import *
import numpy as np
import os

mesh=RectangleMesh(Point(0,0), Point(1, 1), 1, 1, "right/left")

dA = Measure("ds", domain=mesh)
dV = Measure("dx", domain=mesh)


ScalarSpace = FunctionSpace(mesh,'CG',1)
VectorSpace = VectorFunctionSpace(mesh,'CG',1)

test_rho = TestFunction(ScalarSpace)
delta_rho = TrialFunction(ScalarSpace)

v_g_exp = Expression(('1.0','0.0'),degree=2, domain=mesh)
v_g = interpolate(v_g_exp, VectorSpace)


####################  K and D matrix ###########################

k_form = div(delta_rho*v_g)*test_rho*dV
K = assemble(k_form)
D = assemble(k_form)
D.zero()
print K.array()

####################  working  #################################
print D
block = np.ones(9,dtype=np.float_)
print block
rows = np.array([0,1,2],dtype=np.intc)
print rows
cols = np.array([0,1,2],dtype=np.intc)
print cols
D.set(block,rows,cols)
D.apply('add')
print D.array()

####################  not working  ############################
print D
block = np.ones(16,dtype=np.float_)
print block
rows = np.array([0,1,2,3],dtype=np.intc)
print rows
cols = np.array([0,1,2,3],dtype=np.intc)
print cols
D.set(block,rows,cols)
D.apply('add')
print D.array()

###############################################################

(3. How will this react regarding a potential future parallelization?)

Similar question was asked here: https://fenicsproject.org/qa/9017/fem-fct-artificial-diffusion-matrix-construction

asked Aug 16, 2016 by titusf FEniCS Novice (570 points)
edited Aug 16, 2016 by titusf

1 Answer

+6 votes
 
Best answer

The easiest way to do this is using backend functions for accessing matrix data. In particular, you want to avoid looping over matrix entries in the Python code.

The following example shows how you can do this if you have built dolfin with petsc4py.

K = assemble(k_form)

kmat = as_backend_type(K).mat()
dmat = kmat.duplicate()

# copy transpose of kmat into dmat
kmat.transpose(dmat)

# get values from kmat and dmat
i, c, kvals = kmat.getValuesCSR()
_, _, dvals = dmat.getValuesCSR()

# compute values for matrix D
dvals = - reduce(numpy.minimum, (kvals, dvals, 0))

# copy values into matrix D
dmat.setValuesCSR(i, c, dvals)
dmat.assemble()
D = PETScMatrix(dmat)
answered Aug 18, 2016 by Magne Nordaas FEniCS Expert (13,820 points)
selected Jan 7, 2017 by titusf

Helped a lot, thank you.

...