Exact Dirichlet BC

Hallo all,

The code below substitutes the penalized Dirichlet BC with an exact one, while preserving the symmetry of A. With this approach I face the following issue; the introduced zeros increase the nnz of A. Is there a way to consider the new zeros in the sparse format (nnz(Anew)==nnz(Aold))?

border a001(t = 0.0 * pi, 0.5 * pi){ x = cos(t); y = sin(t); label = 0;}
border a002(t = 0.5 * pi, 2.0 * pi){ x = cos(t); y = sin(t); label = 1;}
mesh T = buildmesh(a001(50) + a002(150));

fespace V(T, P1);
varf BND(u, void) = on(1, u = 1.0);
V RHS = 0.0;
RHS[] = BND(0, V, tgv=1.0);

V u, v;
varf lap(u, v) = int2d(T)( dx(u) * dx(v) + dy(u) * dy(v) );
varf rhs(u, v) = int2d(T)( 1.0 * v );
matrix A = lap(V, V);
real[int] b = rhs(0, V);

for (int i = 0; i < V.ndof; i++) {
  if ( abs(RHS[][i] - 1.0) > 1e-10 ) {
    for (int j = 0; j < V.ndof; j++) {
      if ( i != j && abs(RHS[][j] - 1.0) < 1e-10 ) {
        b[i] = b[i] - 50.0 * A(i,j);
      }
    }
  }
}
for (int i = 0; i < V.ndof; i++) {
  if ( abs(RHS[][i] - 1.0) < 1e-10 ) {
    b[i] = 50.0;
    for (int j = 0; j < V.ndof; j++) {
      if ( i != j ) {
        A(i,j) = 0.0;
        A(j,i) = 0.0;
      }
      else {
        A(i,j) = 1.0;
      }
    }
  }
}

u[] = A^-1 * b;
plot(u, fill = true, nbiso=128, wait=true);

V w;
solve poisson(w, v) = int2d(T)( dx(w) * dx(v) + dy(w) * dy(v) ) -  int2d(T)( 1.0 * v ) + on(1, w=50.0);
plot(w, fill = true, nbiso=128, wait=true);

V d = abs(u - w)/abs(u);
plot(d, fill = true, nbiso=128, wait=true);