-- FreeFem++ v4.6 (Thu Apr 2 15:47:38 CEST 2020 - git v4.6) Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // stiffness_matrix.edp 2 : // 3 : // Discussion: 4 : // 5 : // Construct the stiffness matrix for the Poisson equation on the unit disk. 6 : // 7 : // Print the stiffness matrix, and notice that certain entries have unexpected values. 8 : // 9 : // Show that these values are associated with boundary nodes having Dirichlet 10 : // boundary conditions. 11 : // 12 : // Location: 13 : // 14 : // http://people.sc.fsu.edu/~jburkardt/freefem_src/stiffness_matrix/stiffness_matrix.edp 15 : // 16 : // Modified: 17 : // 18 : // 29 June 2015 19 : // 20 : // Author: 21 : // 22 : // John Burkardt 23 : // 24 : cout << "\n"; 25 : cout << "stiffness_matrix:\n"; 26 : cout << " FreeFem++ version:\n"; 27 : cout << "\n"; 28 : cout << " Stiffness matrix definition:\n"; 29 : cout << " K(i,j) = Integral Phi'(i) Phi'(j)\n"; 30 : cout << "\n"; 31 : cout << " For technical reasons, rows of K corresponding ... : to Dirichlet\n"; 32 : cout << " boundary conditions will have a huge number on ... : diagonal.\n"; 33 : cout << "\n"; 34 : cout << " If you don't know this, it may surprise you, or ... : cause your\n"; 35 : cout << " calculations to go astray.\n"; 36 : // 37 : // Define the mesh. 38 : // 39 : mesh Th = square ( 3, 3 ); 40 : // 41 : // Plot the mesh. 42 : // 43 : plot ( Th, ps = "stiffness_matrix_mesh.ps" ); 44 : // 45 : // Define the finite element space. 46 : // 47 : fespace Vh ( Th, P1 ); 48 : Vh uh; 49 : Vh vh; 50 : // 51 : // Define bilinear forms for Dirichlet and Neumann problems. 52 : // 53 : varf bid ( uh, vh ) = int2d ( Th ) ( dx(uh) * dx(vh) + dy(uh) * dy(vh) ) + on ( 1, 2, 3, 4, uh = 0.0 ); 54 : varf bin ( uh, vh ) = int2d ( Th ) ( dx(uh) * dx(vh) + dy(uh) * dy(vh) ); 55 : // 56 : // Create the corresponding stiffness matrices. 57 : // 58 : matrix ad = bid ( Vh, Vh ); 59 : matrix an = bin ( Vh, Vh ); 60 : // 61 : // Print the matrices. 62 : // 63 : cout << "\n"; 64 : cout << " Matrix AD:\n"; 65 : cout << " Stiffness matrix for system with Dirichlet cond ... : itions.\n"; 66 : cout << "\n"; 67 : cout << ad; 68 : cout << "\n"; 69 : cout << " Matrix AN:\n"; 70 : cout << " Stiffness matrix for system with Neumann condit ... : ions.\n"; 71 : cout << "\n"; 72 : cout << an; 73 : // 74 : // Terminate. 75 : // 76 : cout << "\n"; 77 : cout << "stiffness_matrix:\n"; 78 : cout << " Normal end of execution.\n"; 79 : 80 : sizestack + 1024 =1584 ( 560 ) stiffness_matrix: FreeFem++ version: Stiffness matrix definition: K(i,j) = Integral Phi'(i) Phi'(j) For technical reasons, rows of K corresponding to Dirichlet boundary conditions will have a huge number on diagonal. If you don't know this, it may surprise you, or cause your calculations to go astray. -- Square mesh : nb vertices =16 , nb triangles = 18 , nb boundary edges 12 Matrix AD: Stiffness matrix for system with Dirichlet conditions. # HashMatrix Matrix (COO) 0x55d91764f570 # n m nnz half fortran state 16 16 82 0 0 0 0 0 0 1.0000000000000000199e+30 0 1 -0.5 0 5 0 1 0 -0.5 1 1 1.0000000000000000199e+30 1 5 -1 5 0 0 5 1 -1 5 5 4 0 4 -0.5 5 4 -1 4 0 -0.5 4 5 -1 4 4 1.0000000000000000199e+30 1 2 -0.5 1 6 0 2 1 -0.5 2 2 1.0000000000000000199e+30 2 6 -1 6 1 0 6 2 -1 6 6 4 6 5 -1 5 6 -1 2 3 -0.49999999999999988898 2 7 0 3 2 -0.49999999999999988898 3 3 1.0000000000000000199e+30 3 7 -0.5 7 2 0 7 3 -0.5 7 7 1.0000000000000000199e+30 7 6 -0.99999999999999977796 6 7 -0.99999999999999977796 4 9 0 5 9 -1 9 4 0 9 5 -1 9 9 4 4 8 -0.5 9 8 -1 8 4 -0.5 8 9 -1 8 8 1.0000000000000000199e+30 5 10 0 6 10 -1 10 5 0 10 6 -1 10 10 4 10 9 -1 9 10 -1 6 11 0 7 11 -0.5 11 6 0 11 7 -0.5 11 11 1.0000000000000000199e+30 11 10 -0.99999999999999977796 10 11 -0.99999999999999977796 8 13 0 9 13 -0.99999999999999977796 13 8 0 13 9 -0.99999999999999977796 13 13 1.0000000000000000199e+30 8 12 -0.49999999999999988898 13 12 -0.5 12 8 -0.49999999999999988898 12 13 -0.5 12 12 1.0000000000000000199e+30 9 14 0 10 14 -0.99999999999999977796 14 9 0 14 10 -0.99999999999999977796 14 14 1.0000000000000000199e+30 14 13 -0.5 13 14 -0.5 10 15 0 11 15 -0.49999999999999988898 15 10 0 15 11 -0.49999999999999988898 15 15 1.0000000000000000199e+30 15 14 -0.49999999999999988898 14 15 -0.49999999999999988898 Matrix AN: Stiffness matrix for system with Neumann conditions. # HashMatrix Matrix (COO) 0x55d91764fff0 # n m nnz half fortran state 16 16 82 0 0 0 0 0 0 1 0 1 -0.5 0 5 0 1 0 -0.5 1 1 2 1 5 -1 5 0 0 5 1 -1 5 5 4 0 4 -0.5 5 4 -1 4 0 -0.5 4 5 -1 4 4 2 1 2 -0.5 1 6 0 2 1 -0.5 2 2 1.999999999999999778 2 6 -1 6 1 0 6 2 -1 6 6 4 6 5 -1 5 6 -1 2 3 -0.49999999999999988898 2 7 0 3 2 -0.49999999999999988898 3 3 0.99999999999999988898 3 7 -0.5 7 2 0 7 3 -0.5 7 7 1.999999999999999778 7 6 -0.99999999999999977796 6 7 -0.99999999999999977796 4 9 0 5 9 -1 9 4 0 9 5 -1 9 9 4 4 8 -0.5 9 8 -1 8 4 -0.5 8 9 -1 8 8 2 5 10 0 6 10 -1 10 5 0 10 6 -1 10 10 4 10 9 -1 9 10 -1 6 11 0 7 11 -0.5 11 6 0 11 7 -0.5 11 11 1.999999999999999778 11 10 -0.99999999999999977796 10 11 -0.99999999999999977796 8 13 0 9 13 -0.99999999999999977796 13 8 0 13 9 -0.99999999999999977796 13 13 1.999999999999999778 8 12 -0.49999999999999988898 13 12 -0.5 12 8 -0.49999999999999988898 12 13 -0.5 12 12 0.99999999999999988898 9 14 0 10 14 -0.99999999999999977796 14 9 0 14 10 -0.99999999999999977796 14 14 1.999999999999999778 14 13 -0.5 13 14 -0.5 10 15 0 11 15 -0.49999999999999988898 15 10 0 15 11 -0.49999999999999988898 15 15 0.99999999999999977796 15 14 -0.49999999999999988898 14 15 -0.49999999999999988898 stiffness_matrix: Normal end of execution. times: compile 0.004218s, execution 0.000812s, mpirank:0 CodeAlloc : nb ptr 3555, size :473600 mpirank: 0 Ok: Normal End