function []=nbc_ex
%-----------------------------------------------------------------------
%  nbc_ex.m - solves a "generic" elliptic problem in 1D
%  adapted from elliptic_1d.m
%
%  Copyright (c) 2001, Jeff Borggaard, Virginia Tech
%  Version: 1.0
%
%  Usage:    [] = nbc_ex
%
%  Variables:  x_l
%                        Location of the left endpoint
%              x_r
%                        Location of the right endpoint
%              u_l
%                        Dirichlet boundary condition
%              f_1
%                        Flux boundary condition
%
%              elements_mesh
%                        Density of elements in (x_l,x_r)
%
%  Functions:  p_function
%                        function in the second order term ( p>0 )
%              q_function
%                        function in the zeroth order term ( q>=0 )
%-----------------------------------------------------------------------

%-----------------------------------------------------------------------
%  Define "Input" Parameters
%-----------------------------------------------------------------------
  x_l  = 0;
  x_r  = 1;

  u_0  = 1;
  f_1  = 2;

  elements_mesh    = 4;

  variable.element = 'continuous_linear';

  output_type      = 'matlab';

  n_gauss          = 3;  % number of points used in Gaussian integration

%-----------------------------------------------------------------------
%  Geometry Module
%-----------------------------------------------------------------------
  if ( strcmp(variable.element,'continuous_linear') )
    xb(:,1)      = [x_l; x_r];
    e_connb(1,:) = [1 2];
    rho(1)       = elements_mesh;
    [x,e_conn,index_u,index_c] = oned_mesh(xb,e_connb,rho);
  end

  if ( strcmp(variable.element,'continuous_quadratic') )
    xb(:,1)      = [x_l; .5*(x_l+x_r); x_r];
    e_connb(1,:) = [1 2 3];
    rho(1)       = elements_mesh;
    [x,e_conn,index_u,index_c] = oned_mesh(xb,e_connb,rho);
  end

  if ( strcmp(variable.element,'continuous_cubic') )
    xb(:,1)      = [x_l; (x_l+x_r)/3; 2*(x_l+x_r)/3; x_r];
    e_connb(1,:) = [1 2 3 4];
    rho(1)       = elements_mesh;
    [x,e_conn,index_u,index_c] = oned_mesh(xb,e_connb,rho);
  end

  variable.dirichlet_node  = [1];
  variable.dirichlet_value = [u_0];

  [n_elements, nel_dof]    = size(e_conn);
  variable.neumann_element = [n_elements];
  variable.neumann_face    = [2];
  variable.neumann_value   = [p_function(x_r)*f_1];

%-----------------------------------------------------------------------
%  Solver Module
%-----------------------------------------------------------------------
  [n_nodes   , n_dimensions] = size(x     );
  [n_elements, nel_dof     ] = size(e_conn);
  
  n_dirichlet = 0;
  if (isfield(variable,'dirichlet_node'))
    n_dirichlet = n_dirichlet + length(variable.dirichlet_node);
  end
  
  n_neumann = 0;
  if (isfield(variable,'neumann_face'))
    n_neumann = n_neumann + length(variable.neumann_face);
  end
  
  %---------------------------------------------------------------------
  % Determine equation numbers and set up boundary condition information
  %---------------------------------------------------------------------
  n_dir       = 0;
  n_equations = 0;

  for n_nd=1:n_nodes
    is_dir = 0;
    for i=1:length(variable.dirichlet_node)
      if (variable.dirichlet_node(i)==n_nd)
        is_dir = i;
      end
    end  
    % determine whether this degree of freedom corresponds to
    % a Dirichlet boundary condition.
    if (is_dir>0)
      n_dir       = n_dir + 1;
      ide(n_nd)   =-n_dir;
      dir(n_dir)  = variable.dirichlet_value(is_dir);
    else
      n_equations = n_equations + 1;
      ide(n_nd)   = n_equations;
    end
  end

  %---------------------------------------------------------------------
  %  Build the finite element matrices
  %---------------------------------------------------------------------
  A = sparse(n_equations,n_equations);
  b = zeros (n_equations,1          );

  [r,w] = oned_gauss(n_gauss);

  for n_el=1:n_elements
    % compute value of each test function and their spatial derivaties
    % at the integration points
    nodes_local       = e_conn(n_el,:);
    x_local           = x(nodes_local,:);
    [x_g,w_g,phi,p_x] = oned_shape(x_local,r,w);

    % compute the value of functions at the Gauss points
    p_g   = p_function(x_g);
    q_g   = q_function(x_g);
    f_g   = f_function(x_g);

    %--------------------------------------------------------------------
    %  Assemble the weak form of the equations (element contributions)
    %--------------------------------------------------------------------
    A_loc = oned_bilinear(p_g, p_x, p_x, w_g) + ...
            oned_bilinear(q_g, phi, phi, w_g);
    b_loc = oned_f_int   (f_g, phi,      w_g);

    %-----------------------------------------------------------------
    % Assemble contributions into the global system matrix
    %-----------------------------------------------------------------
    for n_t=1:nel_dof
      n_test = ide(nodes_local(n_t));
      if (n_test > 0)  % this is an unknown, fill the row
        for n_u=1:nel_dof
          n_unk = ide(nodes_local(n_u));
          if (n_unk > 0)
            A(n_test,n_unk) = A(n_test,n_unk) + A_loc(n_t,n_u);
          else
            b(n_test) = b(n_test) - A_loc(n_t,n_u)*dir(-n_unk);
          end
        end
        b(n_test) = b(n_test) + b_loc(n_t);
      end
    end
  end

  % add flux contributions
 for n=1:n_neumann
   n_el           = variable.neumann_element(n);
   nodes_el       = e_conn(n_el,:);
   if ( variable.neumann_face(n) == 1 )
     e_local = nodes_el(1);
   else
     e_local = nodes_el(end);
   end
   n_test = ide(e_local);
   b(n_test) = b(n_test) + variable.neumann_value(n);
 end
%  b(end) = b(end) + p_function(x(end,1))*f_1;

  %-----------------------------------------------------------------------
  %  Perform Implicit Solve
  %-----------------------------------------------------------------------
%  figure(5); spy(A);  % view the sparsity pattern in A

  du = A\b;

  %-----------------------------------------------------------------------
  %  Construct the Solution (here, we simply apply the Dirichlet bc's)
  %-----------------------------------------------------------------------
  for n=1:n_nodes
    i = ide(n);
    if (i>0)
      %  get the nodal value out of the linear system solve
      u(n) = du(i);
    else
      %  get the nodal value from the specified Dirichlet bc
      u(n) = dir(-i);
    end
  end

%-----------------------------------------------------------------------
%  Post Processing Module
%-----------------------------------------------------------------------
  %  Write Out Solution
  if ( strcmp(output_type,'matlab') )
    save elliptic1d_output x e_conn u
  elseif ( strcmp(output_type,'none') )
    fprintf('No output is written.')
  else
    fprintf('I omitted these options in 1D.')
  end

  figure(1)
  x_ex = linspace(x(1,1),x(end,1),200);
  u_ex = u_exact(x_ex);

  plot(x,u,':o',x_ex,u_ex)
  legend('Finite Element Solution','Exact Solution')
  xlabel('x'); ylabel('Temperature'); 
  error = norm(u_ex - interp1(x,u,x_ex))
%-----------------------------------------------------------------------
%  Supporting Functions
%-----------------------------------------------------------------------
function f = f_function(x)
  f = 27*x.^4 - 72*x.^3 - 21*x.^2 + 21*x + 6;

function u = u_exact(x)
  u = 3*x.^3 - 4*x.^2 + x + 1;

function p = p_function(x)
  p = (x+1).^2;

function q = q_function(x)
  q = 9*x;

%  oned_bilinear
%  oned_f_int
%  oned_gauss
%  oned_mesh
%  oned_shape