# Exact Dirichlet BC
**URL:** https://community.freefem.org/t/exact-dirichlet-bc/52
**Category:** General Discussion
**Created:** [May 17, 2019, 12:22pm UTC](https://community.freefem.org/t/exact-dirichlet-bc/52 "2019-05-17T12:22:12Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@fotios.kasolis](https://community.freefem.org/u/fotios.kasolis)
#### Post date: [May 17, 2019, 12:22pm UTC](https://community.freefem.org/t/exact-dirichlet-bc/52/1 "2019-05-17T12:22:12Z")
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))`?
```auto
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);
```
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [September 2, 2019, 10:06am UTC](https://community.freefem.org/t/exact-dirichlet-bc/52/2 "2019-09-02T10:06:32Z")
Remark,
this code is too expensive (cost of V.dnof^2)
please done use, and read the doc with tgv
is tgv = -1
a not to bad way to code is:
> 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;}
> int nn=50;
> mesh T = buildmesh(a001(nn) + a002(3\*nn));
>
> fespace V(T, P1);
>
> V u, v , u1;
> varf lap(u, v) = int2d(T)( dx(u) \* dx(v) + dy(u) \* dy(v) ) +on(1,u=50);
> varf rhs(u, v) = int2d(T)( 1.0 \* v ) +on(1,u=50);
> matrix A = lap(V, V,tgv=-1);// remove all row i and put 1 on diag if i is on B.C. (“on” Boundary Condition)
> // C.F. sect 3.3.10 the note in the “on” paragraph.
> // Warning the matrix become A unsymetric
>
> real[int] b = rhs(0, V,tgv=-1); //
>
> // method no symetric
> u = A^-1 \* b;
> // method symetric
> matrix A1 = lap(V, V,tgv=-2);// remove all row and column i and put 1 on diag if i is on B.C.
> // (not in doc to day)
> // the matrix become A1 symetric and Positif definite
>
> varf rhs1(u, v) = int2d(T)( 1.0 \* v ) +on(1,u=0); // with 0 B.C.
> varf rhs2(u, v) = on(1,u=50); // with correct B.C.
> real[int] bc = rhs2(0,V,tgv=-2); // BC. interpolation
> real[int] b1 = rhs1(0,V,tgv=-2); // contribution of RHS
> real[int] b2 = A\*bc; // -BC. contribution but
>
> b1 = b1 -b2; // Warning - BC on BC node due to Id term on B.C
> b1 += 2_bc; // to get corret valeu on B.C. part ( -1+2 = 1 )
> if(verbosity\>2)
> {
> cout \<\< bc \<\< endl;
> cout \<\< A1 \<\< endl;
> cout \<\< b1 \<\< endl;
> }
> u1[] = A1^-1_b1;
>
> plot(u, fill = true, nbiso=128, wait=true,cmm=“u”);
> plot(u1, fill = true, nbiso=128, wait=true,cmm=“u1”);
>
> 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);
---
### Author:  [@fotios.kasolis](https://community.freefem.org/u/fotios.kasolis)
#### Post date: [July 22, 2020, 2:19pm UTC](https://community.freefem.org/t/exact-dirichlet-bc/52/3 "2020-07-22T14:19:56Z")
Dear FreeFEM community,
Although this topic is old - and partially resolved - I think the following should be mentioned. Although my naive implementation is costly, in terms of flops, @frederichecht and @arousta implementations, which are essentially the same, use two sparse matrices and hence, they are inefficient in terms of memory. I guess the way to go here is a direct homogenization of the Dirichlet data. The resulting code should look as simple as
```auto
mesh T = square(10, 10);
fespace V(T, P2);
func uD = y;
// P2 EXTENSION OF DIRICHLET DATA
varf bndDirichlet(u, v) = on(2, 4, u = uD);
V EuD;
EuD[] = bndDirichlet(0, V, tgv = 1.0);
plot(EuD, value = true, fill = true, nbiso=64, wait=true);
// BILINEAR AND LINEAR VARIATIONAL FORMS
varf a(u, v) = int2d(T)( dx(u)*dx(v) + dy(u)*dy(v) ) + on(2, 4, u = 0.0);
varf l(u, v) = int2d(T)( v - dx(EuD)*dx(v) - dy(EuD)*dy(v) ) + on(2, 4, u = 0.0);
// ASSEMBLY
matrix A = a(V, V, tgv = -2); // only one matrix needed
real[int] b = l(0, V); // no need for tgv here, it is multiplied with 0.0
// SOLUTION WITH VANISHING DIRICHLET DATA
V u0;
u0[] = A^-1*b;
plot(u0, value = true, fill = true, nbiso=64, wait=true, cmm = "Solution with vanishing boundary data");
// SOLUTION TO INITIAL PROBLEM
V u = u0 + EuD;
plot(u, value = true, fill = true, nbiso = 64, wait = true, cmm = "Solution to Initial problem");
```
This is actually, how NGSolve implements Dirichlet boundary conditions.
---
### Author:  [@arousta](https://community.freefem.org/u/arousta)
#### Post date: [September 18, 2020, 7:01am UTC](https://community.freefem.org/t/exact-dirichlet-bc/52/4 "2020-09-18T07:01:15Z")
Thanks @fotios.kasolis for sharing the nice idea! I like it 👍