//  stiffness_matrix.edp
//
//  Discussion:
//
//    Construct the stiffness matrix for the Poisson equation on the unit disk.
//
//    Print the stiffness matrix, and notice that certain entries have unexpected values.
//
//    Show that these values are associated with boundary nodes having Dirichlet
//    boundary conditions.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/stiffness_matrix/stiffness_matrix.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 June 2015
//
//  Author:
//
//    John Burkardt
//
cout << "\n";
cout << "stiffness_matrix:\n";
cout << "  FreeFem++ version:\n";
cout << "\n";
cout << "  Stiffness matrix definition:\n";
cout << "    K(i,j) = Integral Phi'(i) Phi'(j)\n";
cout << "\n";
cout << "  For technical reasons, rows of K corresponding to Dirichlet\n";
cout << "  boundary conditions will have a huge number on diagonal.\n";
cout << "\n";
cout << "  If you don't know this, it may surprise you, or cause your\n";
cout << "  calculations to go astray.\n";
//
//  Define the mesh.
//
mesh Th = square ( 3, 3 );
//
//  Plot the mesh.
//
plot ( Th, ps = "stiffness_matrix_mesh.ps" );
//
//  Define the finite element space.
//
fespace Vh ( Th, P1 );
Vh uh;
Vh vh;
//
//  Define bilinear forms for Dirichlet and Neumann problems.
//
varf bid ( uh, vh ) = int2d ( Th ) ( dx(uh) * dx(vh) + dy(uh) * dy(vh) ) 
  + on ( 1, 2, 3, 4, uh = 0.0 );
varf bin ( uh, vh ) = int2d ( Th ) ( dx(uh) * dx(vh) + dy(uh) * dy(vh) );
//
//  Create the corresponding stiffness matrices.
//
matrix ad = bid ( Vh, Vh );
matrix an = bin ( Vh, Vh );
//
//  Print the matrices.
//
cout << "\n";
cout << "  Matrix AD:\n";
cout << "  Stiffness matrix for system with Dirichlet conditions.\n";
cout << "\n";
cout << ad;
cout << "\n";
cout << "  Matrix AN:\n";
cout << "  Stiffness matrix for system with Neumann conditions.\n";
cout << "\n";
cout << an;
//
//  Terminate.
//
cout << "\n";
cout << "stiffness_matrix:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );