Exact Dirichlet BC

General Discussion
1

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);
1 Like
2

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 += 2bc; // 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^-1b1;

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);

1 Like
3

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

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.

2 Likes
4

Thanks @fotios.kasolis for sharing the nice idea! I like it :+1: