#define MAX_ELEM 11744
#define MAX_NODE 2561
#define MAX_K 17933
#define MAX_COMM 1
#define MAX_BOUNDARY 1
#define MAX_U 3280641
#define ELEM_VERTICES 4
#define MESH_DIM 3
#define PARTITIONS 1
#define CELLS 1
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sys/timers.h>
#include "archlib.h"
extern FILE *packfile;
int ARCHself = 0;
char *progname;
struct options options;
int ARCHglobalnode[MAX_NODE];
int ARCHglobalelem[MAX_ELEM];
float ARCHcoord[MAX_NODE][MESH_DIM];
int ARCHvertex[MAX_ELEM][ELEM_VERTICES];
int ARCHcomm[MAX_COMM];
float ARCHcommbuffer[MAX_COMM][3];
int ARCHcommindex[PARTITIONS + 1];
int ARCHmatrixcol[MAX_K + 1];
int ARCHmatrixindex[MAX_NODE + 1];
int ARCHglobalnodes, ARCHmesh_dim;
int ARCHglobalelems, ARCHcorners;
int ARCHsubdomains, ARCHprocessors;
int ARCHelems;
int ARCHnodes, ARCHpriv, ARCHmine;
int ARCHmatrixlen;
int ARCHcholeskylen;
int ARCHcommlen;
int ARCHloopv;
int ARCHloopc1;
int ARCHloopc2;
int ARCHi, ARCHj, ARCHk;
double ARCHiotime = 0.0;
double ARCHcomptime = 0.0;
double ARCHcommtime = 0.0;
double ARCHstarttime;
double ARCHendtime;
struct timespec ARCHtimeval;
struct ip {
  float c_s;
  float c_p;
  float rho;
  float dt;
  float duration;
  float dip;
  float strike;
  float rake;
  float T0;
};
struct properties {
  float p_v;
  float s_v;
  float den;
};
struct source {
  float dip, strike, rake;
  float x0[3];
};
struct source_node {
  int number;
  float x0[3];
  float dist;
};
struct vector {
  double x[3];
};
struct matrix {
  double e[3][3];
};
#define  PI   3.141592654
#define  EPS       0.0001
#define  COEF 0.017453293
#define IDENTITY(x) x
#define XM  20.0
#define YM  20.0
#define ZM -15.0
#define N_obs         2  /* Number of Observation Points */
#define N_obs1    794  /* 1st Node Number */
#define N_obs2    1292  /* 2nd Node Number */
#define N_X_MINUS 192  
#define N_X_PLUS  292
#define N_Y_MINUS 273
#define N_Y_PLUS  293
#define N_Z_MINUS 287
#define N_Z_PLUS  487
#define HX 5.0
#define HY 5.0
#define HZ 5.0
#define Damp_depth -10.0
#define Damp_X_min   4.0
#define Damp_X_max  16.0
#define Damp_Y_min   4.0
#define Damp_Y_max  16.0
struct source s;
struct ip IP;
void KMC_maker();
void prepare();
void iso3d();
void d_shape();
void inv_J();
float phi0(float);
void cal_force(float, int, float *, float *, float *);
void abc_damper();
float area();
void centroid();
int sourcenode[MAX_NODE];
int nodekind[MAX_NODE];
float nodekindf[MAX_NODE];
int forcele[MAX_ELEM];
float elemkindf[3][MAX_ELEM];
float disp[3][MAX_NODE][3], vel[MAX_NODE][3], M[MAX_NODE][3], C[MAX_NODE][3];
float u0[MAX_NODE][3], u2[MAX_NODE][3];
float K[MAX_K][3][3];

void mv3product(nodes, A, Acol, Aindex, v, w)
int nodes;
float A[][3][3];
int *Acol;
int *Aindex;
float v[][3];
float w[][3];
{
  int i;
  int Anext, Alast;
  float vi0, vi1, vi2, sum0, sum1, sum2, value;
  float vcol0, vcol1, vcol2, wcol0, wcol1, wcol2;
  int col;
  
  for (i = 0; i < nodes; i++) {
    w[i][0] = 0.0;
    w[i][1] = 0.0;
    w[i][2] = 0.0;
  }
  for (i = 0; i < nodes; i++) {
    vi0 = v[i][0];
    vi1 = v[i][1];
    vi2 = v[i][2];
    Anext = Aindex[i];
    Alast = Aindex[i + 1];
    sum0 = w[i][0] + A[Anext][0][0] * vi0 + A[Anext][0][1] * vi1 + A[Anext][0][
        2] * vi2;
    sum1 = w[i][1] + A[Anext][1][0] * vi0 + A[Anext][1][1] * vi1 + A[Anext][1][
        2] * vi2;
    sum2 = w[i][2] + A[Anext][2][0] * vi0 + A[Anext][2][1] * vi1 + A[Anext][2][
        2] * vi2;
    Anext++;
    while (Anext < Alast) {
      col = Acol[Anext];
      vcol0 = v[col][0];
      vcol1 = v[col][1];
      vcol2 = v[col][2];
      value = A[Anext][0][0];
      sum0 += value * vcol0;
      wcol0 = w[col][0] + value * vi0;
      value = A[Anext][0][1];
      sum0 += value * vcol1;
      wcol1 = w[col][1] + value * vi0;
      value = A[Anext][0][2];
      sum0 += value * vcol2;
      wcol2 = w[col][2] + value * vi0;
      value = A[Anext][1][0];
      sum1 += value * vcol0;
      wcol0 += value * vi1;
      value = A[Anext][1][1];
      sum1 += value * vcol1;
      wcol1 += value * vi1;
      value = A[Anext][1][2];
      sum1 += value * vcol2;
      wcol2 += value * vi1;
      value = A[Anext][2][0];
      sum2 += value * vcol0;
      w[col][0] = wcol0 + value * vi2;
      value = A[Anext][2][1];
      sum2 += value * vcol1;
      w[col][1] = wcol1 + value * vi2;
      value = A[Anext][2][2];
      sum2 += value * vcol2;
      w[col][2] = wcol2 + value * vi2;
      Anext++;
    }
    w[i][0] = sum0;
    w[i][1] = sum1;
    w[i][2] = sum2;
  }
}


void main(int argc, char **argv)
{
  int i, j, k, iter, timesteps;
  int globalnode;
  int NumNeighbor = 0;
  int NumPlus = 0;
  int NumMinus = 0;
  int globalNumNeighbor = 0;
  int globalNumPlus = 0;
  int disptplus, dispt, disptminus;
  int verticesonbnd;
  int cor[4], bv[3];
  int Np_step;
  float c1[3], c2[3], c3[3], c4[3];
  float time;
  float *Ke[12], Kee[12][12];
  float Me[12], Ce[12], v[12];
  float zeta, const_a, const_b;
  float d_alpha;
  float x[3][4], xc[3];
  float fr;
  struct properties prop;
  int FLOATS_PER_NODE = 7;
  int *ibuf;
  float *fbuf;
  int count;
  int surface_nodes;
  int ARCHloop1;
  int ARCHloop2;
  int ARCHloop3;
  int ARCHloop4;
  int ARCHtemp1;
  
  arch_init(argc, argv, &options);
  if (!options.quiet) fprintf(stderr, "%s: Beginning simulation.\n", argv[0]);
  arch_readnodevector(nodekindf, ARCHnodes);
  arch_readelemvector(elemkindf[0], ARCHelems);
  arch_readelemvector(elemkindf[1], ARCHelems);
  arch_readelemvector(elemkindf[2], ARCHelems);
  IP.dt = 0.04;
  IP.duration = 30.00;
  IP.T0 = 1.0;
  Np_step = 5;
  timesteps = IP.duration / IP.dt + 1;
  s.strike = 0.0;
  s.dip = 90.0;
  s.rake = 0.0;
  s.x0[0] = 9.0;
  s.x0[1] = 9.0;
  s.x0[2] = - 5.0;
  IP.strike = s.strike * PI / 180.0;
  IP.dip = s.dip * PI / 180.0;
  IP.rake = s.rake * PI / 180.0;
  const_a = 0.00533333;
  const_b = 0.06666667;
  zeta = 10.0;
  fr = 0.1;
  for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
    for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
      M[ARCHloopv][ARCHloopc1] = 0.0;
    }
  }
  for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
    for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
      C[ARCHloopv][ARCHloopc1] = 0.0;
    }
  }
  for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
    for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
      disp[0][ARCHloopv][ARCHloopc1] = 0.0;
    }
  }
  for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
    for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
      disp[1][ARCHloopv][ARCHloopc1] = 0.0;
    }
  }
  for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
    for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
      disp[2][ARCHloopv][ARCHloopc1] = 0.0;
    }
  }
  for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
    nodekind[ARCHloopv] = 0;
  }
  for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
    for (ARCHloopc2 = 0; ARCHloopc2 < 3; ARCHloopc2++) {
      for (ARCHloopv = 0; ARCHloopv < ARCHmatrixlen; ARCHloopv++) {
        K[ARCHloopv][ARCHloopc1][ARCHloopc2] = 0.0;
      }
    }
  }
  for (i = 0; i < 12; i++) Ke[i] = Kee[i];
  disptplus = 0;
  dispt = 1;
  disptminus = 2;
  if (ARCHself == 0) {
    fprintf(stderr, "         CASE SUMMARY \n\n");
    fprintf(stderr, " Fault Information\n");
    fprintf(stderr, "   Orientation:  strike: %f\n", s.strike);
    fprintf(stderr, "                    dip: %f\n", s.dip);
    fprintf(stderr, "                   rake: %f\n", s.rake);
    fprintf(stderr, "   Location: (%f, %f, %f), in Km\n", s.x0[0], s.x0[1], s.
        x0[2]);
    fprintf(stderr, "   Ramp rise time: %f, sec\n", IP.T0);
    fprintf(stderr, " Time step: %f, sec\n", IP.dt);
    fprintf(stderr, " Duration: %f, sec\n", IP.duration);
    fflush(stderr);
  }
  fprintf(stdout, "%d\n", N_obs);
  for (i = 0; i < ARCHnodes; i++) {
    c1[0] = ARCHcoord[i][0];
    c1[1] = ARCHcoord[i][1];
    c1[2] = ARCHcoord[i][2];
    if (fabs(c1[0]) < EPS || (fabs(c1[0]) - XM) < EPS) nodekind[i] = 4; else 
        if (fabs(c1[1]) < EPS || (fabs(c1[1]) - YM) < EPS) nodekind[i] = 5;
         else if (fabs(c1[2] - ZM) < EPS) nodekind[i] = 6; else nodekind[i] = 1
        ;
  }
  surface_nodes = 0;
  for (i = 0; i < ARCHmine; i++) {
    c1[0] = ARCHcoord[i][0];
    c1[1] = ARCHcoord[i][1];
    c1[2] = ARCHcoord[i][2];
    if (c1[2] == 0.0) surface_nodes++;
  }
  fprintf(stderr, "For subdomain %d there are %d surface nodes\n", ARCHself, 
      surface_nodes);
  fflush(stderr);
  for (i = 0; i < ARCHelems; i++) {
    for (j = 0; j < 12; j++) {
      Me[j] = 0.0;
      Ce[j] = 0.0;
      v[j] = 0.0;
      for (k = 0; k < 12; k++) Kee[j][k] = 0.0;
    }
    for (ARCHloop1 = 0; ARCHloop1 < ELEM_VERTICES; ARCHloop1++) {
      cor[ARCHloop1] = ARCHvertex[i][ARCHloop1];
    }
    prop.p_v = elemkindf[0][i];
    prop.s_v = elemkindf[1][i];
    prop.den = elemkindf[2][i];
    verticesonbnd = 0;
    for (j = 0; j < 4; j++) if (nodekind[cor[j]] != 1) bv[verticesonbnd++] = j;
    if (verticesonbnd == 3) {
      c1[0] = ARCHcoord[cor[bv[0]]][0];
      c1[1] = ARCHcoord[cor[bv[0]]][1];
      c1[2] = ARCHcoord[cor[bv[0]]][2];
      c2[0] = ARCHcoord[cor[bv[1]]][0];
      c2[1] = ARCHcoord[cor[bv[1]]][1];
      c2[2] = ARCHcoord[cor[bv[1]]][2];
      c3[0] = ARCHcoord[cor[bv[2]]][0];
      c3[1] = ARCHcoord[cor[bv[2]]][1];
      c3[2] = ARCHcoord[cor[bv[2]]][2];
      abc_damper(c1, c2, c3, bv, &prop, Ce);
    }
    for (ARCHloop1 = 0; ARCHloop1 < ELEM_VERTICES; ARCHloop1++) {
      cor[ARCHloop1] = ARCHvertex[i][ARCHloop1];
    }
    c1[0] = ARCHcoord[cor[0]][0];
    c1[1] = ARCHcoord[cor[0]][1];
    c1[2] = ARCHcoord[cor[0]][2];
    c2[0] = ARCHcoord[cor[1]][0];
    c2[1] = ARCHcoord[cor[1]][1];
    c2[2] = ARCHcoord[cor[1]][2];
    c3[0] = ARCHcoord[cor[2]][0];
    c3[1] = ARCHcoord[cor[2]][1];
    c3[2] = ARCHcoord[cor[2]][2];
    c4[0] = ARCHcoord[cor[3]][0];
    c4[1] = ARCHcoord[cor[3]][1];
    c4[2] = ARCHcoord[cor[3]][2];
    KMC_maker(c1, c2, c3, c4, &prop, Ke, Me, Ce);
    for (j = 0; j < 3; j++) {
      x[j][0] = c1[j];
      x[j][1] = c2[j];
      x[j][2] = c3[j];
      x[j][3] = c4[j];
    }
    centroid(x, xc);
    d_alpha = 4.0 * PI * const_a * 0.95 / (prop.s_v + const_b);
    if (xc[0] < Damp_X_min || xc[0] > Damp_X_max) d_alpha = 2.0 * zeta / 100.0 
        * (2.0 * PI * fr);
    if (xc[1] < Damp_Y_min || xc[1] > Damp_Y_max) d_alpha = 2.0 * zeta / 100.0 
        * (2.0 * PI * fr);
    if (xc[2] < Damp_depth) d_alpha = 2.0 * zeta / 100.0 * (2.0 * PI * fr);
    for (j = 0; j < 12; j++) Ce[j] = Ce[j] + d_alpha * Me[j];
    for (ARCHloop1 = 0; ARCHloop1 < ELEM_VERTICES; ARCHloop1++) {
      M[ARCHvertex[i][ARCHloop1]][0] += Me[ARCHloop1 * 3];
      M[ARCHvertex[i][ARCHloop1]][1] += Me[ARCHloop1 * 3 + 1];
      M[ARCHvertex[i][ARCHloop1]][2] += Me[ARCHloop1 * 3 + 2];
    }
    for (ARCHloop1 = 0; ARCHloop1 < ELEM_VERTICES; ARCHloop1++) {
      C[ARCHvertex[i][ARCHloop1]][0] += Ce[ARCHloop1 * 3];
      C[ARCHvertex[i][ARCHloop1]][1] += Ce[ARCHloop1 * 3 + 1];
      C[ARCHvertex[i][ARCHloop1]][2] += Ce[ARCHloop1 * 3 + 2];
    }
    for (ARCHloop1 = 0; ARCHloop1 < ELEM_VERTICES; ARCHloop1++) {
      for (ARCHloop2 = 0; ARCHloop2 < ELEM_VERTICES; ARCHloop2++) {
        if (ARCHvertex[i][ARCHloop1] <= ARCHvertex[i][ARCHloop2]) {
          ARCHtemp1 = ARCHmatrixindex[ARCHvertex[i][ARCHloop1]];
          while (ARCHmatrixcol[ARCHtemp1] != ARCHvertex[i][ARCHloop2]) {
            ARCHtemp1++;
            if (ARCHtemp1 >= ARCHmatrixindex[ARCHvertex[i][ARCHloop2] + 1]) {
              printf("K indexing error!!! %d %d\n", ARCHvertex[i][ARCHloop1], 
                  ARCHvertex[i][ARCHloop2]);
              exit(1);
            }
          }
          for (ARCHloop3 = 0; ARCHloop3 < 3; ARCHloop3++) {
            for (ARCHloop4 = 0; ARCHloop4 < 3; ARCHloop4++) {
              K[ARCHtemp1][ARCHloop3][ARCHloop4] += Kee[ARCHloop1 * 3 + 
                  ARCHloop3][ARCHloop2 * 3 + ARCHloop4];
            }
          }
        }
      }
    }
  }
  ;
  ;
  count = 0;
  for (i = 0; i < ARCHnodes; i++) {
    if (i == ARCHmine) break;
    c1[0] = ARCHcoord[i][0];
    c1[1] = ARCHcoord[i][1];
    c1[2] = ARCHcoord[i][2];
    if (c1[2] == 0.0) {
      nodekind[i] = 3;
      count++;
    }
  }
  for (iter = 1; iter <= timesteps; iter++) {
    time = iter * IP.dt;
    for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
      for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
        u0[ARCHloopv][ARCHloopc1] = 0.0;
      }
    }
    for (i = 0; i < ARCHnodes; i++) {
      globalnode = ARCHglobalnode[i];
      if ((globalnode == N_X_MINUS) || (globalnode == N_X_PLUS) || (globalnode 
          == N_Y_MINUS) || (globalnode == N_Y_PLUS) || (globalnode == N_Z_MINUS
          ) || (globalnode == N_Z_PLUS)) cal_force(time, globalnode, &u0[i][0]
          , &u0[i][1], &u0[i][2]);
    }
    mv3product(ARCHnodes, K, ARCHmatrixcol, ARCHmatrixindex, disp[dispt], disp[
        disptplus]);
    ;
    for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
      for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
        disp[disptplus][ARCHloopv][ARCHloopc1] = IP.dt * IP.dt * (u0[ARCHloopv]
            [ARCHloopc1] - disp[disptplus][ARCHloopv][ARCHloopc1]) + 2 * M[
            ARCHloopv][ARCHloopc1] * disp[dispt][ARCHloopv][ARCHloopc1] + (0.5 
            * IP.dt * C[ARCHloopv][ARCHloopc1] - M[ARCHloopv][ARCHloopc1]) * 
            disp[disptminus][ARCHloopv][ARCHloopc1];
      }
    }
    for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
      for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
        disp[disptplus][ARCHloopv][ARCHloopc1] = disp[disptplus][ARCHloopv][
            ARCHloopc1] / (M[ARCHloopv][ARCHloopc1] + 0.5 * IP.dt * C[ARCHloopv
            ][ARCHloopc1]);
      }
    }
    for (ARCHloopc1 = 0; ARCHloopc1 < 3; ARCHloopc1++) {
      for (ARCHloopv = 0; ARCHloopv < ARCHnodes; ARCHloopv++) {
        vel[ARCHloopv][ARCHloopc1] = 0.5 / IP.dt * (disp[disptplus][ARCHloopv][
            ARCHloopc1] - disp[disptminus][ARCHloopv][ARCHloopc1]);
      }
    }
    if (iter % Np_step == 0) {
      if (ARCHself == 0) fprintf(stderr, "$$$ %d\n", iter);
      for (i = 0; i < ARCHmine; i++) {
        if (ARCHglobalnode[i] == N_obs1) {
          fprintf(stdout, "%14e %14e %14e %14e\n", time, vel[i][0], vel[i][1], 
              vel[i][2]);
          fflush(stdout);
        }
        if (ARCHglobalnode[i] == N_obs2) {
          fprintf(stdout, "%14e %14e %14e %14e\n", time, vel[i][0], vel[i][1], 
              vel[i][2]);
          fflush(stdout);
        }
      }
    }
    i = disptminus;
    disptminus = dispt;
    dispt = disptplus;
    disptplus = i;
  }
  if (!options.quiet) {
    fprintf(stderr, "%s: Done. Terminating the simulation.\n", progname);
  }
  exit(0);
}


void KMC_maker(c1, c2, c3, c4, prop, Ke, Me, Ce)
float *c1, *c2, *c3, *c4, **Ke, *Me, *Ce;
struct properties *prop;
{
  float jicob[3][5], ds[4][4], xyz[8][3], elas[4];
  
  prepare(c1, c2, c3, c4, xyz, elas, prop);
  iso3d(xyz, elas, jicob, ds, Ke, Me);
}


void Enu(e, v, prop)
float *e, *v;
struct properties *prop;
{
  float x;
  
  x = prop->p_v / prop->s_v;
  x = x * x;
  *v = 0.5 * (x - 2) / (x - 1);
  *e = 2 * prop->den * prop->s_v * prop->s_v * (1 + *v);
}


void prepare(c1, c2, c3, c4, xyz, elas, prop)
float *c1, *c2, *c3, *c4;
float xyz[8][3];
float elas[];
struct properties *prop;
{
  int i;
  float e, v;
  
  Enu(&e, &v, prop);
  elas[0] = e / (2.0 * (v + 1.0) * (1.0 - v * 2.0));
  elas[1] = e * v / ((v + 1.0) * (1.0 - v * 2.0));
  elas[2] = e / ((v + 1.0) * 2.0);
  elas[3] = prop->den;
  for (i = 0; i < 3; i++) {
    xyz[0][i] = c1[i];
    xyz[1][i] = c2[i];
    xyz[2][i] = c3[i];
    xyz[3][i] = c4[i];
  }
}


void iso3d(xyz, elas, jicob, ds, Ke, Me)
float xyz[8][3];
float elas[];
float jicob[3][5];
float ds[4][4];
float **Ke;
float *Me;
{
  float det, tt, ts, c1, c2, c3;
  float volume;
  int i, j, m, n, row, column;
  
  d_shape(ds);
  for (i = 0; i < 3; ++i) {
    for (j = 0; j < 3; ++j) {
      det = 0.0;
      for (m = 0; m < 4; ++m) {
        det = det + ds[i][m] * xyz[m][j];
      }
      jicob[j][i] = det;
    }
  }
  inv_J(jicob, &det);
  for (m = 0; m < 4; ++m) {
    for (i = 0; i < 3; ++i) {
      jicob[i][3] = 0.0;
      for (j = 0; j < 3; ++j) {
        jicob[i][3] = jicob[i][3] + jicob[j][i] * ds[j][m];
      }
    }
    for (i = 0; i < 3; ++i) {
      ds[i][m] = jicob[i][3];
    }
  }
  volume = det / 6;
  if (volume <= 0) {
    printf("Warning: volume of the element = %f !\n", volume);
  }
  c1 = elas[0] * volume;
  c2 = elas[1] * volume;
  c3 = elas[2] * volume;
  row = - 1;
  for (m = 0; m < 4; m++) {
    for (i = 0; i < 3; ++i) {
      ++row;
      column = - 1;
      for (n = 0; n <= m; n++) {
        for (j = 0; j < 3; j++) {
          ++column;
          ts = ds[i][m] * ds[j][n];
          if (i == j) {
            ts = ts * c1;
            tt = (ds[0][m] * ds[0][n] + ds[1][m] * ds[1][n] + ds[2][m] * ds[2][
                n]) * c3;
          } else {
            if (m == n) {
              ts = ts * c1;
              tt = 0;
            } else {
              ts = ts * c2;
              tt = ds[j][m] * ds[i][n] * c3;
            }
          }
          Ke[row][column] = Ke[row][column] + ts + tt;
        }
      }
    }
  }
  tt = elas[3] * volume / 4.0;
  for (i = 0; i < 12; i++) Me[i] = tt;
  for (i = 0; i < 12; ++i) for (j = 0; j <= i; j++) Ke[j][i] = Ke[i][j];
}


void d_shape(ds)
float ds[4][4];
{
  ds[0][0] = - 1;
  ds[1][0] = - 1;
  ds[2][0] = - 1;
  ds[0][1] = 1;
  ds[1][1] = 0;
  ds[2][1] = 0;
  ds[0][2] = 0;
  ds[1][2] = 1;
  ds[2][2] = 0;
  ds[0][3] = 0;
  ds[1][3] = 0;
  ds[2][3] = 1;
}


void inv_J(a, det)
float a[3][5];
float *det;
{
  float d1;
  float c[3][3];
  int i, j;
  
  c[0][0] = a[1][1] * a[2][2] - a[2][1] * a[1][2];
  c[0][1] = a[0][2] * a[2][1] - a[0][1] * a[2][2];
  c[0][2] = a[0][1] * a[1][2] - a[0][2] * a[1][1];
  c[1][0] = a[1][2] * a[2][0] - a[1][0] * a[2][2];
  c[1][1] = a[0][0] * a[2][2] - a[0][2] * a[2][0];
  c[1][2] = a[0][2] * a[1][0] - a[0][0] * a[1][2];
  c[2][0] = a[1][0] * a[2][1] - a[1][1] * a[2][0];
  c[2][1] = a[0][1] * a[2][0] - a[0][0] * a[2][1];
  c[2][2] = a[0][0] * a[1][1] - a[0][1] * a[1][0];
  *det = a[0][0] * c[0][0] + a[0][1] * c[1][0] + a[0][2] * c[2][0];
  d1 = 1.0 / *det;
  for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) a[i][j] = c[i][j] * d1;
}


float cal_area(c1, c2, c3)
float *c1, *c2, *c3;
{
  float a, b, c;
  float x2, y2, z2;
  float p;
  float area;
  
  x2 = (c1[0] - c2[0]) * (c1[0] - c2[0]);
  y2 = (c1[1] - c2[1]) * (c1[1] - c2[1]);
  z2 = (c1[2] - c2[2]) * (c1[2] - c2[2]);
  a = sqrt(x2 + y2 + z2);
  x2 = (c3[0] - c2[0]) * (c3[0] - c2[0]);
  y2 = (c3[1] - c2[1]) * (c3[1] - c2[1]);
  z2 = (c3[2] - c2[2]) * (c3[2] - c2[2]);
  b = sqrt(x2 + y2 + z2);
  x2 = (c1[0] - c3[0]) * (c1[0] - c3[0]);
  y2 = (c1[1] - c3[1]) * (c1[1] - c3[1]);
  z2 = (c1[2] - c3[2]) * (c1[2] - c3[2]);
  c = sqrt(x2 + y2 + z2);
  p = (a + b + c) / 2;
  area = sqrt(p * (p - a) * (p - b) * (p - c));
  if (area <= 0) {
    fprintf(stderr, "Error: negative area or zero area %f\n", p);
    exit(- 1);
  }
  return area;
}


int area_dir(c1, c2, c3)
float *c1, *c2, *c3;
{
  if (fabs(c1[0] - c2[0]) <= EPS && fabs(c2[0] - c3[0]) <= EPS) return 0; else 
      if (fabs(c1[1] - c2[1]) <= EPS && fabs(c2[1] - c3[1]) <= EPS) return 1;
       else if (fabs(c1[2] - c2[2]) <= EPS && fabs(c2[2] - c3[2]) <= EPS) 
      return 2;
}


void abc_damper(c1, c2, c3, bv, prop, Ce)
int bv[3];
float *c1, *c2, *c3;
float Ce[12];
struct properties *prop;
{
  int i, m0, ad;
  float area;
  
  area = cal_area(c1, c2, c3);
  ad = area_dir(c1, c2, c3);
  if (ad == 0) for (i = 0; i < 3; i++) {
    m0 = 3 * bv[i];
    Ce[m0] = Ce[m0] + prop->p_v * prop->den * area / 3;
    Ce[m0 + 1] = Ce[m0 + 1] + prop->s_v * prop->den * area / 3;
    Ce[m0 + 2] = Ce[m0 + 2] + prop->s_v * prop->den * area / 3;
  } else if (ad == 1) for (i = 0; i < 3; i++) {
    m0 = 3 * bv[i];
    Ce[m0] = Ce[m0] + prop->s_v * prop->den * area / 3;
    Ce[m0 + 1] = Ce[m0 + 1] + prop->p_v * prop->den * area / 3;
    Ce[m0 + 2] = Ce[m0 + 2] + prop->s_v * prop->den * area / 3;
  } else if (ad == 2) for (i = 0; i < 3; i++) {
    m0 = 3 * bv[i];
    Ce[m0] = Ce[m0] + prop->s_v * prop->den * area / 3;
    Ce[m0 + 1] = Ce[m0 + 1] + prop->s_v * prop->den * area / 3;
    Ce[m0 + 2] = Ce[m0 + 2] + prop->p_v * prop->den * area / 3;
  }
}


float phi0(float t)
{
  float value;
  
  if (t <= IP.T0) {
    value = 0.5 / PI * (2.0 * PI * t / IP.T0 - sin(2.0 * PI * t / IP.T0));
    return value;
  } else return 1.0;
}


void centroid(float x[3][4], float x0[3])
{
  int i;
  
  for (i = 0; i < 3; i++) x0[i] = (x[i][0] + x[i][1] + x[i][2] + x[i][3]) / 4;
}

struct matrix smdt();

void cal_force(float t, int node, float *force1, float *force2, float *force3)
{
  float f0;
  float M0, h0;
  struct matrix m;
  
  m = smdt();
  M0 = 2.25e+3;
  h0 = 5.0;
  f0 = M0 / (2 * h0) * phi0(t);
  *force1 = 0;
  *force2 = 0;
  *force3 = 0;
  if (node == N_X_MINUS) {
    *force1 = - m.e[0][0] * f0 * h0 / HX;
    *force2 = - m.e[0][1] * f0 * h0 / HX;
    *force3 = - m.e[0][2] * f0 * h0 / HX;
  } else if (node == N_X_PLUS) {
    *force1 = m.e[0][0] * f0 * h0 / HX;
    *force2 = m.e[0][1] * f0 * h0 / HX;
    *force3 = m.e[0][2] * f0 * h0 / HX;
  } else if (node == N_Y_MINUS) {
    *force1 = - m.e[1][0] * f0 * h0 / HY;
    *force2 = - m.e[1][1] * f0 * h0 / HY;
    *force3 = - m.e[1][2] * f0 * h0 / HY;
  } else if (node == N_Y_PLUS) {
    *force1 = m.e[1][0] * f0 * h0 / HY;
    *force2 = m.e[1][1] * f0 * h0 / HY;
    *force3 = m.e[1][2] * f0 * h0 / HY;
  } else if (node == N_Z_MINUS) {
    *force1 = - m.e[2][0] * f0 * h0 / HZ;
    *force2 = - m.e[2][1] * f0 * h0 / HZ;
    *force3 = - m.e[2][2] * f0 * h0 / HZ;
  } else if (node == N_Z_PLUS) {
    *force1 = m.e[2][0] * f0 * h0 / HZ;
    *force2 = m.e[2][1] * f0 * h0 / HZ;
    *force3 = m.e[2][2] * f0 * h0 / HZ;
  }
}


struct matrix smdt()
{
  int i, j;
  struct vector n, t;
  struct matrix m;
  
  t.x[0] = cos(IP.rake) * sin(IP.strike) - sin(IP.rake) * cos(IP.strike) * cos(
      IP.dip);
  t.x[1] = cos(IP.rake) * cos(IP.strike) + sin(IP.rake) * sin(IP.strike) * cos(
      IP.dip);
  t.x[2] = sin(IP.rake) * sin(IP.dip);
  n.x[0] = cos(IP.strike) * sin(IP.dip);
  n.x[1] = - sin(IP.strike) * sin(IP.dip);
  n.x[2] = cos(IP.dip);
  for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) m.e[i][j] = n.x[i] * t.x[j] +
       n.x[j] * t.x[i];
  return m;
}