[scripts] properly clear everything outside the PoU support

parent b704fc69
......@@ -281,13 +281,12 @@ class PouLocalizedCoerciveReductor(DistributedReductorBase):
if (II, JJ) in lhs_global:
local_op = lhs_global[II, JJ]
mat = local_op.assemble(mu=mu).matrix.copy()
# This should not matter mathematically, but worsens the estimate quite a bit.
# clear_rows_and_cols(
# mat,
# row_indices=np.where(pou_and_grads_are_zero[II])[0],
# unit_diag=False,
# prune=True,
# col_indices=np.where(pou_and_grads_are_zero[JJ])[0])
clear_rows_and_cols(
mat,
row_indices=np.where(pou_and_grads_are_zero[II])[0],
unit_diag=False,
prune=True,
col_indices=[])
lhs_pwn[ii, jj] = NumpyMatrixOperator(mat, local_op.source.id, local_op.range.id)
lhs_pwn = BlockOperator(lhs_pwn)
rhs_global = self.fom.operators['rhs']
......@@ -295,6 +294,7 @@ class PouLocalizedCoerciveReductor(DistributedReductorBase):
for ii, II in enumerate(patch):
local_rhs = rhs_global[II]
mat = local_rhs.assemble(mu=mu).matrix.copy()
mat[0, np.where(pou_and_grads_are_zero[II])[0]] *= 0
rhs_patch[0, ii] = NumpyMatrixOperator(mat, local_rhs.source.id, local_rhs.range.id)
rhs_patch = BlockOperator(rhs_patch)
# build patch product
......@@ -307,7 +307,7 @@ class PouLocalizedCoerciveReductor(DistributedReductorBase):
prod_patch = bmat(prod_patch)
prod_patch.eliminate_zeros()
# restrict patch product to support of PoU
prod_DoFs_to_clear = np.where(np.hstack([pou_is_zero[ss] for ss in patch]))[0]
prod_DoFs_to_clear = np.where(np.hstack([pou_and_grads_are_zero[ss] for ss in patch]))[0]
prod_patch = clear_rows_and_cols(prod_patch, prod_DoFs_to_clear, unit_diag=True, prune=True)
prod_patch = NumpyMatrixOperator(prod_patch)
# reconstruct solution
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment