### Apply boundary condition to a manipulated matrix

200
views
1
6 months ago by
Hallo,

I want to assemble matrices, then I want to do some algebraic manipulations (with petsc4py), and finally I want to solve the problem with Dirichlet boundary conditions. However, I cannot apply the boundary conditions (Error: Unable to successfully call PETSc function 'VecSetValuesLocal').

from dolfin import *
import numpy as np
from petsc4py import PETSc

##### DATA #####
mesh=UnitSquareMesh(4,4)
f = Constant(1.0)
u_D = Constant(0.0)

##### FUNCTION SPACES #####
RT_elem = FiniteElement('RT', mesh.ufl_cell(),1)
S_elem = FiniteElement('P', mesh.ufl_cell(),1)
Y_elem = FiniteElement('DG', mesh.ufl_cell(),1)
X = FunctionSpace(mesh,MixedElement([S_elem,RT_elem]))
Y = FunctionSpace(mesh,Y_elem)

##### FORMS #####
y1 = TrialFunction(Y)
y2 = TestFunction(Y)
u,p = TrialFunctions(X)
rhs = (f*y2)*dx
dummy = inner(Constant(1),y2)*dx

##### MATRICES AND VECTORS #####
G = PETScMatrix()
B = PETScMatrix()
F = PETScVector()
assemble_system(a, rhs, A_tensor=G, b_tensor=F)
assemble_system(b, dummy, A_tensor=B)

G_mat = as_backend_type(G).mat()
B_mat = as_backend_type(B).mat()
F_vec = as_backend_type(F).vec()

##### COMPUTE MATRIX K = B^T * G * B #####
K = B_mat.transposeMatMult(G_mat*B_mat)
K = PETScMatrix(K)

##### COMPUTE VECTOR RHS_vec = B^T G^T F #####
RHS_vec = PETSc.Vec().createWithArray(np.empty(B.size(1)))
temp_vec = F_vec.copy()
G_mat.multTranspose(F_vec,temp_vec)
B_mat.multTranspose(temp_vec,RHS_vec)
RHS_vec = PETScVector(RHS_vec)

##### SOLVE K x = RHS_vec with BC #####
bc = DirichletBC(X.sub(0), u_D, 'on_boundary')
bc.apply(K,RHS_vec)

x = Function(X)
xVec = x.vector()
solve(K,xVec,RHS_vec)​

Where does the error come from and is there a possibility to circumvent the problem?

Best regards,
Johannes

Community: FEniCS Project

3
6 months ago by
The error seems to come from the fact that the objects K and RHS_vec do not have the local-to-global mappings needed for the BC application.  A possible work-around might be to use linear algebra objects that are created from function spaces.  The following code at least executes without error, although it feels a bit kludgy to assemble a made-up form for K (and you might want to check the results for correctness...):

from dolfin import *
import numpy as np
from petsc4py import PETSc

##### DATA #####
mesh=UnitSquareMesh(4,4)
f = Constant(1.0)
u_D = Constant(0.0)

##### FUNCTION SPACES #####
RT_elem = FiniteElement('RT', mesh.ufl_cell(),1)
S_elem = FiniteElement('P', mesh.ufl_cell(),1)
Y_elem = FiniteElement('DG', mesh.ufl_cell(),1)
X = FunctionSpace(mesh,MixedElement([S_elem,RT_elem]))
Y = FunctionSpace(mesh,Y_elem)

##### FORMS #####
y1 = TrialFunction(Y)
y2 = TestFunction(Y)
u,p = TrialFunctions(X)
rhs = (f*y2)*dx
dummy = inner(Constant(1),y2)*dx

##### MATRICES AND VECTORS #####
G = PETScMatrix()
B = PETScMatrix()
F = PETScVector()
assemble_system(a, rhs, A_tensor=G, b_tensor=F)
assemble_system(b, dummy, A_tensor=B)

G_mat = as_backend_type(G).mat()
B_mat = as_backend_type(B).mat()
F_vec = as_backend_type(F).vec()

##### COMPUTE MATRIX K = B^T * G * B #####

####################### MODIFIED ##############################################
#K = B_mat.transposeMatMult(G_mat*B_mat)
#K = PETScMatrix(K)
asdf = TrialFunction(X)
jkl = TestFunction(X)
K = PETScMatrix()
assemble(inner(asdf,jkl)*dx, tensor=K)
B_mat.transposeMatMult(G_mat*B_mat, result=as_backend_type(K).mat())
###############################################################################

##### COMPUTE VECTOR RHS_vec = B^T G^T F #####

####################### MODIFIED ##############################################
#RHS_vec = PETSc.Vec().createWithArray(np.empty(B.size(1)))
RHS_vec = as_backend_type(Function(X).vector()).vec()
###############################################################################

temp_vec = F_vec.copy()
G_mat.multTranspose(F_vec,temp_vec)
B_mat.multTranspose(temp_vec,RHS_vec)
RHS_vec = PETScVector(RHS_vec)

##### SOLVE K x = RHS_vec with BC #####
bc = DirichletBC(X.sub(0), u_D, 'on_boundary')
bc.apply(K,RHS_vec)

x = Function(X)
xVec = x.vector()
solve(K,xVec,RHS_vec)
​
Great, this workaround works perfect. Thanks a lot.
written 6 months ago by Johannes