/* fem1d.c -- 1次元 Poisson 方程式を有限要素法で解く
 *
 *   curl -O https://m-katsurada.sakura.ne.jp/program/fem/fem1d.c
 *   cc -o fem1d fem1d.c
 *   ./fem1d
 *   gnuplot
 *   gnuplot> plot "fem1d.out" with lp
 *   gnuplot> quit
 */

/*
 * 【境界値問題】
 * Ω=(a,b) は数直線上の領域 (開区間)で、
 * その境界 Γ は Γ1, Γ2 と二つの部分からなる。
 * Ω 上与えられた関数 f に対して、
 *
 *              -u''=f  in Ω
 *
 *              u=0  on Γ1,  ∂u/∂n=0  on Γ2
 *
 * を満たす u=u(x) を求めよ。ここで n は Γ の外向き単位法線ベクトル。
 *
 * f≡1,Γ1={0},Γ2={1} つまり -u''=1, u(0)=u'(1)=0 とする
 *  u(x)=x(2-x)/2=-x^2/2+x
 *
 * 出力 fem1d.out を gnuplot で処理すると解のグラフが描ける。
 *  gnuplot> plot "fem1d.out" with linespoints
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef double *vector;
typedef vector *matrix;

double  f(double);
vector  new_vector(int);
matrix  new_matrix(int, int); 
void    assem(int, int, int, matrix, vector, vector x, int ibc[]);
void    ecm(int, vector, double [2][2], double [2]);
void    output(int, vector);
void    output2(int nnode, vector x, vector u);
void    solve(matrix, vector, int);

/* 方程式の非同次項 */
double f(double x)
{
    return 1.0;
}

/*      main program
 *      the finite element method for the Poisson equation
 */

int main(void)
{
    /*
     * nnode        節点の総数 (m)
     * nelmt        有限要素の総数
     * nbc          基本境界条件を課す節点の総数
     * x[]          節点の座標
     * ibc[]        基本境界条件を課す節点の番号
     * am[][]       全体係数行列
     * fm[]         全体ベクトル
     *
     * 注意: 節点番号は 0 から、要素番号は1からつける。よって
     *   節点番号 0,1,...,nnode-1
     *   要素番号 1,...,nelmt
     * のような範囲にわたる
     *
     */

    int ie, nnode, nelmt, nbc;
    double h;
    vector x, fm;
    matrix am;
    int ibc[2];

    /* 総要素数 */
    nelmt = 100;
    /* 節点の総数 */
    nnode = nelmt + 1;
    /* ディリクレ境界にある節点の個数 */
    nbc = 1;
    /* メモリーの確保 */
    x = new_vector(nnode);
    fm = new_vector(nnode);
    am = new_matrix(nnode, nnode);
    /* 節点の座標 (ここでは等分割する) */
    h = 1.0 / nelmt;
    for (ie = 0; ie < nnode; ie++)
      x[ie] = ie * h;
    /* ディリクレ境界にある節点の番号 */
    ibc[0] = 0;

    assem(nnode, nelmt, nbc, am, fm, x, ibc);
    solve(am, fm, nnode);
    output(nnode, fm);
    output2(nnode, x, fm);
    return 0;
}

/* "直接剛性法" --- 連立1次方程式を組み上げる */
void assem(int nnode, int nelmt, int nbc, matrix am, vector fm,
           vector x, int ibc[2])
{
    int i, j, ie, ii, jj;
    /* the direct stiffness method (直接剛性法) */
    double ae[2][2], fe[2];
    /* initial clear */
    for (i = 0; i < nnode; i++) {
        fm[i] = 0.0;
        for (j = 0; j < nnode; j++)
            am[i][j] = 0.0;
    }
    /* assemblage of total matrix and vector (全体行列,全体ベクトルの組立) */
    for (ie = 1; ie <= nelmt; ie++) {
        ecm(ie, x, ae, fe);
        for (i = 0; i < 2; i++) {
            ii = ie + i - 1;
            fm[ii] += fe[i];
            for (j = 0; j < 2; j++) {
                jj = ie + j - 1;
                am[ii][jj] += ae[i][j];
            }
        }
    }
    /* the homogeneous Dirichlet condition (同次Dirichlet境界条件) */
    /* 番号 ibc[0]..ibc[nbc-1] の節点での値は既知であるから、
     * その番号(=:ii)の行と列は対角成分以外 := 0、対角成分 := 1 とする。
     */
    for (i = 0; i < nbc; i++) {
        ii = ibc[i];
        fm[ii] = 0.0;
        for (j = 0; j < nnode; j++)
            am[ii][j] = am[j][ii] = 0.0;
        am[ii][ii] = 1.0;
    }
}

/* 要素係数行列, 要素自由ベクトルを計算する */
void ecm(int k, vector x, double ae[2][2], double fe[2])
{
  double f0 = f(x[k-1]), f1 = f(x[k]);
  /* k番目の要素 [x_{k-1},x_{k}] の長さ */
  double d = x[k] - x[k-1];

  /* 要素係数行列 A_e=(<Li,Lj>)
   *
   *   これは講義でもやったように
   *
   *             1        [ 1  -1]
   *  A_k = ------------- |      |
   *         x[k]-x[k-1]  [-1   1]
   */

  ae[0][0] =   1 / d; ae[0][1] = - 1 / d;
  ae[1][0] = - 1 / d; ae[1][1] =   1 / d;

  /* 要素自由項ベクトル
   *
   * f_k=((f,L0),(f,L1))^T
   *
   * 準備として
   *  (L0,L0)=(L1,L1)=(x[k]-x[k-1])∫    x^2 dx=(x[k]-x[k-1])/3
   *                                [0,1]
   *  (L0,L1)=(L1,L0)=(x[k]-x[k-1])∫    x(1-x) dx=(x[k]-x[k-1])/6
   *                                [0,1]
   * f を 1 次補間する
   *  f≒f0 L0+f1 L1
   *
   * 1次補間に基づいた近似
   *  (f,L0)≒f0(L0,L0)+f1(L1,L0)=(x[k]-x[k-1])(2f0+f1)/6
   *  (f,L1)≒f0(L0,L1)+f1(L1,L1)=(x[k]-x[k-1])(f0+2f1)/6
   *
   */

  fe[0] = d * (2 * f0 + f1) / 6;
  fe[1] = d * (f0 + 2 * f1) / 6;
}

/* Gauss の消去法により連立一次方程式を解く */
void solve(matrix a, vector f, int m)
{
  int m1, i, j, k;
  double aa;
  /* forward elimination; */
  m1 = m - 1;
  for (i = 0; i < m1; i++) {
    for (j = i + 1; j < m; j++) {
      aa = a[j][i] / a[i][i];
      f[j] -= aa * f[i];
      for (k = i + 1; k < m; k++)
        a[j][k] -= aa * a[i][k];
    }
  }
  /* backward substitution */
  f[m-1] /= a[m-1][m-1];
  for (i = m - 2; i >= 0; i--) {
    for (j = i + 1; j < m; j++)
      f[i] -= a[i][j] * f[j];
    f[i] /= a[i][i];
  }
}

void output(int nnode, vector fm)
{
    int i;
    /* output of approximate nodal values of u */
    printf("nodal values of u (節点での u の値)\n");
    for (i = 0; i < 3; i++)
        printf("   i      u    ");
    for (i = 0; i < nnode; i++) {
        if (i % 3 == 0) printf("\n");
        printf("%4d %11.3e", i, fm[i]);
    }
    printf("\n");
}

void output2(int nnode, vector x, vector u)
{
  int i;
  FILE *fp;
  for (i = 0; i < nnode; i++)
    printf("%f %f\n", x[i], u[i]);
  if ((fp = fopen("fem1d.out", "w")) != NULL) {
    for (i = 0; i < nnode; i++)
      fprintf(fp, "%f %f\n", x[i], u[i]);
    fclose(fp);
  }
}

vector new_vector(int n)
{
    return malloc(sizeof(double) * n);
}

void del_vector(vector a)
{
    free(a);
}

matrix new_matrix(int m, int n)
{
  int i;
  matrix a;
  if ((a = malloc(m * sizeof(vector))) == NULL) {
    return NULL;
  }
  for (i = 0; i < m; i++) {
    if ((a[i] = new_vector(n)) == NULL) {
      while (--i >= 0) free(a[i]);
      free(a);
      return NULL;
    }
  }
  return a;
}