/*\
|*|
|*|     ******   ******   *       *   *    *
|*|     *        *        *       *    *  *
|*|     *****    *****    *       *     **
|*|     *        *        *       *    *  *
|*|     *        ******   *****   *   *    *
|*|
|*|
|*|     (c) 1994 - 1996 Andreas Kolms
|*|
|*|     Die theoretischen Grundlagen zu diesem rudimentaeren Finite-Elemente-
|*|     Programm werden im Rahmen der Vorlesung "Numerische Methoden in der
|*|     Statik" unter "1. Einfuehrung in die Methode der finiten Elemente"
|*|     behandelt. Sehr hilfreich bei der Erstellung war das Buch "Finite-
|*|     Elemente-Methoden" von Klaus-Juergen Bathe, berichtigter Nachdruck,
|*|     Springer-Verlag, Berlin, Heidelberg 1990. Die vorliegende Version ist
|*|     in ANSI C programmiert ohne Ausnutzung von Spracherweiterungen (siehe
|*|     z. B. "Programmieren in C" von Kernighan/Ritchie, 2. Ausgabe, Hanser
|*|     Verlag, Muenchen, Wien 1990). Es ist auch eine FORTRAN 77-Version
|*|     verfuegbar.
|*|
|*|     "Felix" ist nur fuer Lehrzwecke gedacht. Es kann frei kopiert,
|*|     veraendert und erweitert werden. Jeder Kopie bitte eine Originalversion
|*|     des Quellcodes, d. h. diese Datei, beifuegen.
|*|
|*|     Das Finite-Elemente-Modell wird ueber eine Eingabedatei definiert, die
|*|     Berechnungsergebnisse werden in eine Ausgabedatei geschrieben und die
|*|     Kommunikation erfolgt mittels "stdin", "stdout" und "stderr".
\*/

#define VERSION "Felix 0.4"
#define LINE 81

/*
 *  Standard Definitionsdateien
 */
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <stdarg.h>
#include <float.h>

/*
 *  Funktionsprototypen
 */
void arguments(int argc, char *argv[], char **ifile, char **ofile);
void line(FILE *in, char *mess, char *string, int max);
void read(FILE *in, char *mess, char *format, void *arg);
void print(FILE *out, char *format, ...);
void error(char *format, ...);
char *cmem(int size);
int *imem(int size);
double *dmem(int size);
int **iimem(int size1, int size2);
double **ddmem(int size1, int size2);
void nodes(int m23, int nn, int *neq, int ***id, double ***x, FILE *in,
    FILE *out);
void loads(int m23, int nn, int neq, int **id, double **r, FILE *in, FILE *out);
void elements(int m23, int nn, int **id, int ne, int neq, int **ntype,
    int ***nele, double ***ele, int **mht, FILE *in, FILE *out);
void eltype(int m23, int ntype, int *nne, int *npe);
void colht(int m23, int nne, int *nele, int **id, int neq, int *mht);
void addresses(int neq, int *m, int *nwk, FILE *out);
void assembly(int m23, int **id, double **x, int ne, int *ntype, int **nele,
    double **ele, int *mda, int nwk, const double eps, double **a);
void truss(int m23, double *x1, double *x2, double e, double a, int ie,
    const double eps, double *s);
void quad(double *x1, double *x2, double *x3, double *x4, int iq, int nint,
    double thic, double e, double pr, int ie, const double eps, double *s);
void gauss(int nint, double *xg, double *wgt);
void matm(double e, double pr, int iq, const double eps, double *c);
void stdm(double *x1, double *x2, double *x3, double *x4, double r, double s,
    int ie, int ic, const double eps, double d[][8], double *det, double *rad);
void addele(int m23, int nne, int *nele, int **id, double *s, int *mda,
    double *a);
void colsol(int ktr, int neq, int *mda, const double eps, double *a, double *v);
void displacements(int m23, int nn, int neq, int **id, double *v, FILE *out);

/*\
|*|     Hauptfunktion
|*|
|*|     Input:
|*|
|*|        argc       = Anzahl der Argumente
|*|        argv[argc] = Argumente
|*|
|*|     Return:
|*|
|*|        0 = normale Beendigung
|*|        1 = Fehler
\*/
main(int argc, char *argv[])
{
    const double eps = sqrt(DBL_EPSILON);
    char head[LINE], *ifile, *ofile;
    int m23, nn, neq, **id, ne, *ntype, **nele, *m, nwk, i;
    double **x, *r, **ele, *a;
    FILE *in, *out;

    arguments(argc, argv, &ifile, &ofile);
    if ((in = fopen(ifile, "r")) == NULL) error("Eingabedatei \"%s\"", ifile);
    if ((out = fopen(ofile, "w")) == NULL) error("Ausgabedatei \"%s\"", ofile);
    print(out, VERSION"\n");

    line(in, "Kopfzeile", head, LINE);
    read(in, "Modell", "%d", &m23);
    read(in, "Anzahl der Knoten", "%d", &nn);
    print(out, "\n%s\n\n%dD Modell\n\nAnzahl der Knoten\n%d\n",head, m23, nn);
    if (m23 != 2 && m23 != 3) error("m23 != 2 && m23 != 3");
    if (nn <= 0) error("nn <= 0");
/*
 *  Knotendaten lesen und Freiheitsgrade bestimmen
 */
    nodes(m23, nn, &neq, &id, &x, in, out);
/*
 *  Einzellasten lesen und Lastvektor aufstellen
 */
    loads(m23, nn, neq, id, &r, in, out);

    read(in, "Anzahl der Elemente", "%d", &ne);
    print(out, "\nAnzahl der Elemente\n%d\n", ne);
    if (ne < 1) error("ne < 1");
/*
 *  Elementdaten lesen und die Hoehen der aktiven Spalten der
 *  Gesamtsteifigkeitsmatrix bestimmen
 */
    elements(m23, nn, id, ne, neq, &ntype, &nele, &ele, &m, in, out);
/*
 *  Adressen der Diagonalelemente der reduzierten Gesamtsteifigkeitsmatrix aus
 *  den Spaltenhoehen berechnen
 */
    addresses(neq, m, &nwk, out);
/*
 *  Elementsteifigkeitsmatrizen bestimmen und in die reduzierte
 *  Gesamtsteifigkeitsmatrix einbauen
 */
    assembly(m23, id, x, ne, ntype, nele, ele, m, nwk, eps, &a);
/*
 *  Gleichungssystem loesen
 */
    colsol(0, neq, m, eps, a, r);
/*
 *  Verschiebungen ausgeben
 */
    displacements(m23, nn, neq, id, r, out);
/*
 *  Speicher freigeben und beenden
 */
    for (i = 0; i < nn; ++i) {free(id[i]); free(x[i]);}
    for (i = 0; i < ne; ++i) {free(nele[i]); free(ele[i]);}
    free(id); free(x); free(r); free(ntype); free(nele); free(ele); free(m); 
    free(a);
    print(out, "\nProgramm beendet\n");
    return 0;
}

/*\
|*|     Namen fuer Eingabe- und Ausgabedatei aus den Argumenten fischen
|*|
|*|     Input:
|*|
|*|         argc       = Anzahl der Argumente
|*|         argv[argc] = Argumente aus dem Programmaufruf
|*|
|*|     Output:
|*|
|*|         ifile      = Name der Eingabedatei
|*|         ofile      = Name der Ausgabedatei
\*/
void arguments(int argc, char *argv[], char **ifile, char **ofile)
{
    if (!argc || argc == 1) {
        *ifile = cmem(LINE);
        print(stdout, "Name der Eingabedatei:\n");
        line(stdin, "Name der Eingabedatei", *ifile, LINE);
        *ofile = cmem(LINE);
        print(stdout, "Name der Ausgabedatei:\n");
        line(stdin, "Name der Ausgabedatei ", *ofile, LINE);
    } else if (argc != 3) {
        error("falsche Anzahl Argumente (argc = %d)", argc);
    } else {
        *ifile = argv[1];
        *ofile = argv[2];
    }
}

/*\
|*|     Zeile einlesen
|*|
|*|     Input:
|*|
|*|         in          = Eingabestrom
|*|         message     = Text fuer Fehlermeldung
|*|         string[max] = Zeichenvektor
|*|         max         = Groesse des Zeichenvektors
|*|
|*|     Output:
|*|
|*|         arg         = Zeichenvektor
\*/
void line(FILE *in, char *mess, char *string, int max)
{
    int i, c;

    --max;
    for (i=0; i<max; ++i) {
        c = fgetc(in);
        if (c == EOF)
            error("%s nicht lesbar", mess);
        else 
            break;

        string[i] = c;
    }
    string[i] = '\0';
}

/*\
|*|     Formatierte Eingabe
|*|
|*|     Input:
|*|
|*|         in      = Eingabestrom
|*|         message = Text fuer Fehlermeldung
|*|         format  = Format-Zeichenkette (wie fscanf)
|*|
|*|     Output:
|*|
|*|         arg     = Argument
\*/
void read(FILE *in, char *mess, char *format, void *arg)
{
    if (fscanf(in, format, arg) == EOF) error("%s nicht lesbar", mess);
}

/*\
|*|     Formatierte Ausgabe
|*|
|*|     Input:
|*|
|*|         out    = Ausgabestrom
|*|         format = Format-Zeichenkette (wie fprintf)
|*|         ...    = variable Argumentliste
\*/
void print(FILE *out, char *format, ...)
{
    va_list ags;

    va_start(ags, format);
    vfprintf(out, format, ags);
    va_end(ags);
}

/*\
|*|     Programm mit Fehlermeldung und -status beenden
|*|
|*|     Input:
|*|
|*|         format = Format-Zeichenkette (wie fprintf)
|*|         ...    = variable Argumentliste
\*/
void error(char *format, ...)
{
    va_list ags;

    fprintf(stderr, "\n\n** STOP: ");
    va_start(ags, format);
    vfprintf(stderr, format, ags);
    va_end(ags);
    fprintf(stderr, "\n\a");
    exit(1);
}

/*\
|*|     Speicher fuer "char"-Feld anfordern
|*|
|*|     Input:
|*|
|*|         size    = Laenge des Feldes
|*|
|*|     Return:
|*|
|*|         s[size] = Zeiger auf den Speicherbereich
\*/
char *cmem(int size)
{
    char *s;

    if (size < 1)
        error("size = %d", size);
    else if ((s = (char *) malloc(size * sizeof(char))) == NULL)
        error("kein Speicher (size = %d)", size);
    return s;
}

/*\
|*|     Speicher fuer "int"-Feld anfordern
|*|
|*|     Input:
|*|
|*|         size    = Laenge des Feldes
|*|
|*|     Return:
|*|
|*|         s[size] = Zeiger auf den Speicherbereich
\*/
int *imem(int size)
{
    int *s;

    if (size < 1)
        error("size = %d", size);
    else if ((s = (int *) malloc(size * sizeof(int))) == NULL)
        error("kein Speicher (size = %d)", size);
    return s;
}

/*\
|*|     Speicher fuer "double"-Feld anfordern
|*|
|*|     Input:
|*|
|*|         size    = Laenge des Feldes
|*|
|*|     Return:
|*|
|*|         s[size] = Zeiger auf den Speicherbereich
\*/
double *dmem(int size)
{
    double *s;

    if (size < 1)
        error("size = %d", size);
    else if ((s = (double *) malloc(size * sizeof(double))) == NULL)
        error("kein Speicher (size = %d)", size);
    return s;
}

/*\
|*|     Speicher fuer zweidimensionales "int"-Feld anfordern
|*|
|*|     Input:
|*|
|*|         size1           = Laenge des Feldes
|*|         size2           = Breite des Feldes
|*|
|*|     Return:
|*|
|*|         s[size1][size2] = Zeiger auf den Speicherbereich
\*/
int **iimem(int size1, int size2)
{
    int i, **s;

    if (size1 < 1)
        error("size1 = %d", size1);
    else if (size2 < 1)
        error("size2 = %d", size2);
    else if ((s = (int **) malloc(size1 * sizeof(int *))) == NULL)
        error("kein Speicher (size1 = %d)", size1);
    for (i = 0; i < size1; ++i)
        if ((s[i] = (int *) malloc(size2 * sizeof(int))) == NULL)
            error("kein Speicher (size2 = %d)", size2);
    return s;
}

/*\
|*|     Speicher fuer zweidimensionales "double"-Feld anfordern
|*|
|*|     Input:
|*|
|*|         size1           = Laenge des Feldes
|*|         size2           = Breite des Feldes
|*|
|*|     Return:
|*|
|*|         s[size1][size2] = Zeiger auf den Speicherbereich
\*/
double **ddmem(int size1, int size2)
{
    int i;
    double **s;

    if (size1 < 1)
        error("size1 = %d", size1);
    else if (size2 < 1)
        error("size2 = %d", size2);
    else if ((s = (double **) malloc(size1 * sizeof(double *))) == NULL)
        error("kein Speicher (size1 = %d)", size1);
    for (i = 0; i < size1; ++i)
        if ((s[i] = (double *) malloc(size2 * sizeof(double))) == NULL)
            error("kein Speicher (size2 = %d)", size2);
    return s;
}

/*\
|*|     Knotendaten lesen und Freiheitsgrade bestimmen
|*|
|*|     Input:
|*|
|*|         m23         = 2 = zweidimensionales Modell
|*|                     = 3 = dreidimensionales Modell
|*|         nn          = Anzahl der Knoten
|*|         in          = Eingabestrom
|*|         out         = Ausgabestrom
|*|
|*|     Output:
|*|
|*|         neq         = Anzahl der Freiheitsgrade
|*|         id[nn][m23] = Freiheitsgradnummern
|*|         x[nn][m23]  = Knotenkoordinaten
\*/
void nodes(int m23, int nn, int *neq, int ***id, double ***x, FILE *in,
    FILE *out)
{
    int i, j, n;

    print(out, "\nKnoten Randbedingungen Koordinaten\n");
    *id = iimem(nn, m23);
    *x = ddmem(nn, m23);
    for (j = 0; j < nn; ++j) {
        read(in, "Knotennummer", "%d", &n);
        print(out, "%d  ", n);
        if (j != (n - 1)) error("Unordnung bei den Knotendaten");
        for (i = 0; i < m23; ++i) {
            read(in, "Randbedingungen", "%d", &(*id)[j][i]);
            print(out, " %d", (*id)[j][i]);            /* 0 = fest, !0 = frei */
        }
        print(out, "  ");
        for (i = 0; i < m23; ++i) {
            read(in, "Knotenkoordinaten", "%lf", &(*x)[j][i]);
            print(out, " %g", (*x)[j][i]);
        }
        print(out, "\n");
    }

    *neq = 0;                                     /* Freiheitsgrade bestimmen */
    for (j = 0; j < nn; ++j)
        for (i = 0; i < m23; ++i)
            if ((*id)[j][i]) (*id)[j][i] = ++(*neq);   /* neuer Freiheitsgrad */
            else             (*id)[j][i] = 0;           /* kein Freiheitsgrad */

    print(out, "\nKnoten Freiheitsgrade\n");
    for (j = 0; j < nn; ++j) {
        print(out, "%d  ", j + 1);
        for (i = 0; i < m23; ++i) print(out, " %d", (*id)[j][i]);
        print(out, "\n");
    }
    if (!(*neq)) error("keine Freiheitsgrade");
    print(out, "\nAnzahl der Freiheitsgrade\n%d\n", *neq);
}

/*\
|*|     Einzellasten lesen und Lastvektor aufstellen
|*|
|*|     Input:
|*|
|*|         m23         = 2 = zweidimensionales Modell
|*|                     = 3 = dreidimensionales Modell
|*|         nn          = Anzahl der Knoten
|*|         new         = Anzahl der Freiheitsgrade
|*|         id[nn][m23] = Freiheitsgradnummern
|*|         in          = Eingabestrom
|*|         out         = Ausgabestrom
|*|
|*|     Output:
|*|
|*|         r[neq]      = Lastvektor
\*/
void loads(int m23, int nn, int neq, int **id, double **r, FILE *in, FILE *out)
{
    int i, nl, ln, li;
    double fl;

    read(in, "Anzahl der Einzellasten", "%d", &nl);
    print(out, "\nAnzahl der Einzellasten\n%d\n", nl);
    if (nl < 1) error("nl < 1");

    *r = dmem(neq);
    for (i = 0; i < neq; ++i) (*r)[i] = 0.0;
    print(out, "\nKnoten Richtung Last\n");
    for (i = 0; i < nl; ++i) {
        read(in, "Knoten", "%d", &ln);
        read(in, "Richtung", "%d", &li);
        read(in, "Last", "%lf", &fl);
        print(out, "%d %d %g\n", ln, li, fl);
        if (li > m23 || li < 1 || ln < 1 || ln > nn)
            error("Fehler bei den Einzellasten");
        if (id[--ln][--li] > 0 && id[ln][li] <= neq) (*r)[id[ln][li] - 1] += fl;
    }
}

/*\
|*|     Elementdaten lesen und die Hoehen der aktiven Spalten der
|*|     Gesamtsteifigkeitsmatrix bestimmen
|*|
|*|     Input:
|*|
|*|         m23           = 2 = zweidimensionales Modell
|*|                       = 3 = dreidimensionales Modell
|*|         nn            = Anzahl der Knoten
|*|         id[nn][m23]   = Freiheitsgradnummern
|*|         ne            = Anzahl der Elemente
|*|         neq           = Anzahl der Freiheitsgrade
|*|         in            = Eingabestrom
|*|         out           = Ausgabestrom
|*|
|*|     Output:
|*|
|*|         ntype[ne]     = Elementtypen
|*|         nele[ne][nne] = Elementdefinitionen
|*|         ele[ne][npe]  = Elementparameter
|*|         mht[neq]      = Hoehen der aktiven Spalten (ohne Diagonale)
\*/
void elements(int m23, int nn, int **id, int ne, int neq, int **ntype,
    int ***nele, double ***ele, int **mht, FILE *in, FILE *out)
{
    int i, j, n, nne, npe;
    const int nne_max = 4, npe_max = 3;

    *ntype = imem(ne);
    *nele = iimem(ne, nne_max);
    *ele = ddmem(ne, npe_max);
    *mht = imem(neq + 1);                               /* (--> mda[neq + 1]) */

    for (i = 0; i < neq; ++i) (*mht)[i] =(int) 0.0;                 /* Anfangswert */

    for (i = 0; i < ne; ++i) {
        read(in, "Elementtyp", "%d", &(*ntype)[i]);

        eltype(m23, (*ntype)[i], &nne, &npe);  /* Knoten- und Parameteranzahl */
        if (nne > nne_max) error("nne_max zu klein");
        if (npe > npe_max) error("npe_max zu klein");

        read(in, "n", "%d", &n);
        for (j = 0; j < nne; ++j)
            read(in, "Elementdefinition", "%d", &(*nele)[i][j]);
        if (i != (n - 1))
            error("Unordnung bei der Elementdefinition (Element %d)", n);
        for (j = 0; j < nne; ++j)
            if ((*nele)[i][j] < 1 || (*nele)[i][j] > nn)
                error("Fehler bei der Elementdefinition (Element %d)", i + 1);

                                                      /* Hoehen aktualisieren */
        colht(m23, nne, (*nele)[i], id, neq, *mht);
        for (j = 0; j < npe; ++j)
            read(in, "Elementparameter", "%lf", &(*ele)[i][j]);
        print(out, "\nElement Elementtyp Elementdefinition Elementparameter\n");
        print(out, "%d %d  ", n, (*ntype)[i]);
        for (j = 0; j < nne; ++j) print(out, " %d", (*nele)[i][j]);
        print(out, "  ");
        for (j = 0; j < npe; ++j) print(out, " %g", (*ele)[i][j]);
        print(out, "\n");
    }
}

/*\
|*|     Knoten- und Parameteranzahl dem Elementtyp zuordnen
|*|
|*|     Input:
|*|
|*|         m23   = 2 = zweidimensionales Modell
|*|               = 3 = dreidimensionales Modell
|*|         ntype = Elementtyp
|*|
|*|     Output:
|*|
|*|         nne   = Anzahl der Knoten des Elements
|*|         npe   = Anzahl der Parameter des Elements (oder NULL)
\*/
void eltype(int m23, int ntype, int *nne, int *npe)
{
    if (ntype == 1) {                                          /* Stabelement */
        *nne = 2;
        if (npe != NULL) *npe = 2;

    } else if (ntype == 2 || ntype == 3 || ntype == 4) {    /* Viereckelement */
        if (m23 != 2) error("Elementtyp %d nur als"
            " zweidimensionales Modell moeglich", ntype);
        *nne = 4;
        if (npe != NULL) *npe = 3;
    } else                                             /* Unbekanntes Element */
        error("unbekannter Elementtyp (ntype = %d", ntype);
}

/*\
|*|     Hoehen der aktiven Spalten der Gesamtsteifigkeitsmatrix aktualisieren
|*|
|*|     Input:
|*|
|*|         m23         = 2 = zweidimensionales Modell
|*|                     = 3 = dreidimensionales Modell
|*|         nne         = Anzahl der Knoten des Elements
|*|         nele[nne]   = Knotennummern des Elements
|*|         id[nn][m23] = Freiheitsgradnummern
|*|         neq         = Anzahl der Freiheitsgrade
|*|         mht[neq]    = Hoehen der aktiven Spalten (ohne Diagonale)
|*|
|*|     Output:
|*|
|*|         mht[neq]    = aktualisierte Hoehen der aktiven Spalten
\*/
void colht(int m23, int nne, int *nele, int **id, int neq, int *mht)
{
    int i, j, idmin, me, ls, idi;

    ls = neq;                       /* Niedrigster Freiheitsgrad des Elements */
    for (i = 0; i < nne; ++i) {
        for (j = 0; j < m23; ++j) {
            idmin = id[nele[i] - 1][j];
            if (!idmin ) continue;
            if (idmin < 0 || idmin > neq) error("idmin < 0 || idmin > neq");
            if (idmin < ls) ls = idmin;
        }
    }

    for (i = 0; i < nne; ++i) {                /* Spaltenhoehen aktualisieren */
        for (j = 0; j < m23; ++j) {
            if (!(idi = id[nele[i] - 1][j])) continue;
            me = idi - ls;
            if (me > mht[--idi]) mht[idi] = me;
        }
    }                                         /* mht[0] ist immer gleich Null */
}

/*\
|*|     Adressen der Diagonalelemente der reduzierten Gesamtsteifigkeitsmatrix
|*|     aus den Spaltenhoehen berechnen
|*|
|*|     Input:
|*|
|*|         neq      = Anzahl der Gleichungen
|*|         m[neq]   = Hoehen der aktiven Spalten (ohne Diagonale)
|*|         out      = Ausgabestrom
|*|
|*|     Output:
|*|
|*|         m[neq+1] = Adressen der Diagonalelemente
|*|         nwk      = Anzahl der Elemente innerhalb der Kontur der Matrix
\*/
void addresses(int neq, int *m, int *nwk, FILE *out)
{
    int i, mm, mht;

    m[0] = 0;
    mht = m[1];
    m[1] = 1;
    for (i = 2; i <= neq; ++i) {
        mm = mht;
        mht = m[i];
        m[i] = m[i - 1] + mm + 1;
    }
    *nwk = m[neq] - m[0];
    print(out, "\nAnzahl der Elemente der reduzierten"
        " Gesamtsteifigkeitsmatrix\n%d\n", *nwk);
}

/*
|*|     Elementsteifigkeitsmatrizen bestimmen und in die reduzierte
|*|     Gesamtsteifigkeitsmatrix einbauen
|*|
|*|     Input:
|*|
|*|         m23           = 2 = zweidimensionales Modell
|*|                       = 3 = dreidimensionales Modell
|*|         id[nn][m23]   = Freiheitsgradnummern
|*|         x[nn][m23]    = Knotenkoordinaten
|*|         ne            = Anzahl der Elemente
|*|         ntype[ne]     = Elementtypen
|*|         nele[ne][nne] = Elementdefinitionen
|*|         ele[ne][npe]  = Elementparameter
|*|         mda[neq+1]    = Adressen der Diagonalelemente
|*|         nwk           = Anzahl der Elemente innerhalb der Kontur der Matrix
|*|         eps           = Genauigkeit
|*|
|*|     Output:
|*|
|*|         a(nwk)     = reduzierte Gesamtsteifigkeitsmatrix
|*|
|*|     Es wird die Huellentechnik (skyline storage) verwendet, d. h.
|*|     Nullelemente ausserhalb der Kontur der Gesamtsteifigkeitsmatrix werden
|*|     nicht gespeichert.
\*/
void assembly(int m23, int **id, double **x, int ne, int *ntype, int **nele,
    double **ele, int *mda, int nwk, const double eps, double **a)
{
    int i, nne;
#define s_max 36
    double s[s_max];

    *a = dmem(nwk);
    for (i = 0; i < nwk; ++i) (*a)[i] = 0.0;                   /* Anfangswert */

    for (i = 0; i < ne; ++i) {

        eltype(m23, ntype[i], &nne, NULL);                    /* Knotenanzahl */
        if (m23*nne*(m23*nne+1)/2 > s_max) error("s_max zu klein");

        if (ntype[i] == 1)                       /* Elementsteifigkeitsmatrix */
            truss(m23, x[nele[i][0] - 1], x[nele[i][1] - 1], ele[i][0],
                ele[i][1], i + 1, eps, s);

        else if (ntype[i] == 2 || ntype[i] == 3 || ntype[i] == 4)
            quad(x[nele[i][0] - 1], x[nele[i][1] - 1], x[nele[i][2] - 1],
                x[nele[i][3] - 1], ntype[i] - 1, 2, ele[i][0], ele[i][1],
                ele[i][2], i, eps, s);

                                        /* Elementsteifigkeitsmatrix einbauen */
        addele(m23, nne, nele[i], id, s, mda, *a);
    }
}

/*\
|*|     Stabelement (Fachwerkstab, d. h. uebertraegt nur axiale Kraefte)
|*|
|*|     Input:
|*|
|*|         m23      = 2 = zweidimensionales Modell
|*|                  = 3 = dreidimensionales Modell
|*|         x1[m23]  = Knotenkoordinaten des ersten Knotens
|*|         x2[m23]  = Knotenkoordinaten des zweiten Knotens
|*|         e        = Elastizitaetsmodul
|*|         a        = Querschnittsflaeche
|*|         ie       = Elementnummer
|*|         eps      = Genauigkeit
|*|
|*|     Output:
|*|
|*|         s[10/21] = Elementsteifigkeitsmatrix
|*|
|*|     o--------o
|*|     1        2
|*|
|*|         | s[0] s[1] s[2]  s[3]  s[4]  s[5]  |
|*|         |      s[6] s[7]  s[8]  s[9]  s[10] |
|*|         |           s[11] s[12] s[13] s[14] |   |  S -S |
|*|     s = |                 s[15] s[16] s[17] | = |       | ea/l**3 ,
|*|         |                       s[18] s[19] |   | -S  S |
|*|         |                             s[20] |
|*|
|*|         |   (x2-x1)**2   (x2-x1)(y2-y1) (x2-x1)(z2-z1) |
|*|     S = | (x2-x1)(y2-y1)   (y2-y1)**2   (y2-y1)(z2-z1) |
|*|         | (x2-x1)(z2-z1) (y2-y1)(z2-z1)   (z2-z1)**2   |
|*|
|*|     bzw.
|*|
|*|         | s[0] s[1] s[2] s[3] |
|*|         |      s[4] s[5] s[6] |   |  S -S |
|*|     s = |           s[7] s[8] | = |       | ea/l**3 ,
|*|         |                s[9] |   | -S  S |
|*|
|*|         |   (x2-x1)**2   (x2-x1)(y2-y1) |
|*|     S = | (x2-x1)(y2-y1)   (y2-y1)**2   |
\*/
void truss(int m23, double *x1, double *x2, double e, double a, int ie,
    const double eps, double *s)
{
    int i, j, k;
    double  xl2, xl, d[3], st[6], xx, yy;

    xl2 = 0.0;
    for (i = 0; i < m23; ++i) {
        d[i] = x1[i] - x2[i];
        xl2 += d[i] * d[i];
    }
    if (xl2 < eps) error("Stab hat die Laenge Null (Element %d)", ie);
    xl = sqrt(xl2);
    xx = e * a * xl;
    for (i = 0; i < m23; ++i) {
        st[i] = d[i] / xl2;
        st[i + m23] = -st[i];
    }
    k = 0;
    for (i = 0; i < (2 * m23); ++i) {
        yy = st[i] * xx;
        for (j = i; j < (2 * m23); ++j) s[k++] = st[j] * yy;
    }
}


/*\ 
|*|  Balkenelement mit 3 Freiheitsgraden pro Knoten
|*|
|*|
|*|
|*|
|*|
|*|
|*|
\*/
void beam(int m23, double *x1, double *x2, double e, double a, int ie,  
          const double eps, double *s)
{
	int i, j, k;
	double x23=0.0;
}

/*\
|*|     Isoparametrisches Viereckelement
|*|
|*|     Input:
|*|
|*|         x1[2] = Knotenkoordinaten des ersten Knotens
|*|         x2[2] = Knotenkoordinaten des zweiten Knotens
|*|         x3[2] = Knotenkoordinaten des dritten Knotens
|*|         x4[2] = Knotenkoordinaten des vierten Knotens
|*|         iq    = 1 = ebener Spannungszustand
|*|               = 2 = ebener Verzerrungszustand
|*|               = 3 = Rotationssymmetrie
|*|         nint  = Integrationsordnung (1 ... 4)
|*|         thic  = Dicke (bzw. Winkel im Bogenmass fuer iq = 3)
|*|         e     = Elastizitaetsmodul
|*|         pr    = Querkontraktionszahl
|*|         ie    = Elementnummer
|*|         eps   = Genauigkeit
|*|
|*|     Output:
|*|
|*|         s[36] = Elementsteifigkeitsmatrix
|*|
|*|              /   T           /1   /1   T
|*|         s =  |  d c d  dB =  |    |   d c d  thic det(J)  dr ds =
|*|             /B              /-1  /-1
|*|
|*|           | s[0] s[1] s[2]  s[3]  s[4]  s[5]  s[6]  s[7]  |
|*|           |      s[8] s[9]  s[10] s[11] s[12] s[13] s[14] |
|*|           |           s[15] s[16] s[17] s[18] s[19] s[20] |
|*|           |                 s[21] s[22] s[23] s[24] s[25] |
|*|         = |                       s[26] s[27] s[28] s[29] |
|*|           |                             s[30] s[31] s[32] |
|*|           |                                   s[33] s[34] |
|*|           |                                         s[35] |
|*|
|*|                       (-1, 1)  2          1  ( 1, 1)
|*|         y                       o--------o
|*|         |                       |        |
|*|         |                       |        |
|*|         |                       |        |
|*|         |                       o--------o
|*|         |             (-1,-1)  3          4  ( 1,-1)
|*|         +--- x
\*/
void quad(double *x1, double *x2, double *x3, double *x4, int iq, int nint,
    double thic, double e, double pr, int ie, const double eps, double *s)
{
    int i, j, k, l, ic, il, jk, kl, lx, ly;
    double c[10], xg[4], wgt[4], d[4][8], det, rad = 1.0, wt, db[4], stiff;

    if (iq != 1 && iq != 2 && iq != 3) error("falscher Aufruf von \"quad\"");
    else if (iq == 3) ic = 4;
    else ic = 3;

    matm(e, pr, iq, eps, c);                  /* Spannungs-Verzerrungs-Matrix */

    gauss(nint, xg, wgt);               /* Stuetzstellen und Gewichtsfaktoren */

    for (i = 0; i < 36; ++i) s[i] = 0.0;         /* Elementsteifigkeitsmatrix */
    for (lx = 0; lx < nint; ++lx) {
        for (ly = 0; ly < nint; ++ly) {
                    /* Verzerrungs-Verschiebungs-Matrix, Determinante, Radius */
            stdm(x1, x2, x3, x4, xg[lx], xg[ly], ie, ic, eps, d, &det, &rad);
            wt = wgt[lx] * wgt[ly] * thic * rad * det;
            il = 0;
            for (j = 0; j < 8; ++j) {
                jk = -1;
                for (k = 0; k < ic; ++k) {
                    jk += k + 1;
                    kl = jk;
                    db[k] = 0.0;
                    for (l = 0; l < k; ++l)
                        db[k] += d[l][j] * c[kl--];
                    for (l = k; l < ic; ++l) {
                        db[k] += d[l][j] * c[kl];
                        kl += l + 2;
                    }
                }
                for (i = j; i < 8; ++i) {
                    stiff = 0.0;
                    for (l = 0; l < ic; ++l) stiff += db[l] * d[l][i];
                    s[il++] += stiff * wt;
                }
            }
        }
    }
}

/*\
|*|     Konstanten fuer die numerische Intergration (Gausssche Quadratur)
|*|
|*|     Input:
|*|
|*|         nint      = Integrationsordnung (1 ... 4)
|*|
|*|     Output:
|*|
|*|         xg[nint]  = Stuetzstellen
|*|         wgt[nint] = Gewichtsfaktoren
\*/
void gauss(int nint, double *xg, double *wgt)
{
    if (nint == 1)
        xg[0] =  0.0, wgt[0] =  2.0;
    else if (nint == 2)
        xg[0] = -0.5773502691896, wgt[0] =  1.0,
        xg[1] =  0.5773502691896, wgt[1] =  1.0;
    else if (nint == 3)
        xg[0] = -0.7745966692415, wgt[0] = 0.5555555555556,
        xg[1] =  0.0,             wgt[1] = 0.8888888888889,
        xg[2] =  0.7745966692415, wgt[2] = 0.5555555555556;
    else if (nint == 4)
        xg[0] = -0.8611363115941, wgt[0] = 0.3478548451375,
        xg[1] = -0.3399810435849, wgt[1] = 0.6521451548625,
        xg[2] =  0.3399810435849, wgt[2] = 0.6521451548625,
        xg[3] =  0.8611363115941, wgt[3] = 0.3478548451375;
    else
        error("falscher Aufruf von \"gauss\"");
}

/*\
|*|     Spannungs-Verzerrungs-Matrix (Stoffmatrix, Elastizitaetsmatrix) fuer
|*|     Scheiben und rotationssymmetrische Systeme bestimmen
|*|
|*|     Input:
|*|
|*|         e     = Elastizitaetsmodul
|*|         pr    = Querkontraktionszahl
|*|         iq    = 1 = ebener Spannungszustand (nc = 6)
|*|               = 2 = ebener Verzerrungszustand (nc = 6)
|*|               = 3 = Rotationssymmetrie (nc = 10)
|*|         eps   = Genauigkeit
|*|
|*|     Output:
|*|
|*|         c[nc] = Elastizitaetsmatrix (Hookesches Gesetz, isotrop)
|*|
|*|         | c[0] c[2] c[5] c[9] |
|*|         |      c[1] c[4] c[8] |
|*|     c = |           c[3] c[7] |
|*|         |                c[6] |
\*/
void matm(double e, double pr, int iq, const double eps, double *c)
{
    double c1, c2, c3;

    if (e < eps)
        error("Elastizitaetsmodul %g nicht zulaessig)", e);
    else if (pr < 0.0 || pr > ((iq == 1) ? 0.5 : 0.5 - eps))
        error("Querkontraktionszahl %g nicht zulaessig)", pr);

    if (iq == 1) {                                 /* Ebener Spannungszustand */
        c1 = e / (1.0 - pr * pr);
        c2 = pr * c1;
        c3 = (1.0 - pr) * c1;

    } else {              /* Ebener Verzerrungszustand und Rotationssymmetrie */
        c3 = e / (pr + 1.0);
        c2 = c3 * pr / (1.0 - pr * 2.0);
        c1 = c2 + c3;
    }
    c[0] = c1;
    c[1] = c1;
    c[2] = c2;
    c[3] = c3 * 0.5;
    c[4] = 0.0;
    c[5] = 0.0;

    if (iq == 3) {                                      /* Rotationssymmetrie */
        c[6] = c1;
        c[7] = 0.0;
        c[8] = c2;
        c[9] = c2;
    }
}

/*\
|*|     Verzerrungs-Verschiebungs-Matrix, Determinante der Jacobi-Matrix und
|*|     gegebenenfalls Radius fuer Viereckelemente im Punkt r, s der
|*|     Parameterebene mittels bilinearer Ansatzfunktionen bestimmen
|*|
|*|     Input:
|*|
|*|         x1[2]    = Knotenkoordinaten des ersten Knotens
|*|         x2[2]    = Knotenkoordinaten des zweiten Knotens
|*|         x3[2]    = Knotenkoordinaten des dritten Knotens
|*|         x4[2]    = Knotenkoordinaten des vierten Knotens
|*|         r        = natuerliche Koordinate (-1 <= r <= 1)
|*|         s        = natuerliche Koordinate (-1 <= s <= 1)
|*|         ie       = Elementnummer
|*|         ic       = 3 = Scheibenelement
|*|                  = 4 = rotationssymmetrisches Element
|*|         eps      = Genauigkeit
|*|
|*|     Output:
|*|
|*|         d[ic][8] = Verzerrungs-Verschiebungs-Matrix
|*|         det      = Determinante der Jacobi-Matrix
|*|         rad      = Radius (nur fuer ic = 4)
\*/
void stdm(double *x1, double *x2, double *x3, double *x4, double r, double s,
    int ie, int ic, const double eps, double d[][8], double *det, double *rad)
{
    int i, j, k, k1, k2;
    double rp, sp, rm, sm, h[4], p[2][4], xj[2][2], xjj;

    if (ic != 3 && ic != 4) error("falscher Aufruf von \"stdm\"");

    rp = 1.0 + r;
    sp = 1.0 + s;
    rm = 1.0 - r;
    sm = 1.0 - s;

    p[0][0] =  0.25 * sp;                           /* Ableitungen nach r     */
    p[0][1] = -p[0][0];                             /* der Ansatzfunktionen   */
    p[0][2] = -0.25 * sm;                           /* fuer r und s auswerten */
    p[0][3] = -p[0][2];

    p[1][0] =  0.25 * rp;                           /* Ableitungen nach s     */
    p[1][1] =  0.25 * rm;                           /* der Ansatzfunktionen   */
    p[1][2] = -p[1][1];                             /* fuer r und s auswerten */
    p[1][3] = -p[1][0];

    for (i = 0; i < 2; ++i)                                  /* Jacobi-Matrix */
        for (j = 0; j < 2; ++j)
            xj[i][j] = p[i][0] * x1[j] +
                       p[i][1] * x2[j] +
                       p[i][2] * x3[j] +
                       p[i][3] * x4[j];

    *det = xj[0][0] * xj[1][1] -            /* Determinante der Jacobi-Matrix */
           xj[0][1] * xj[1][0];
    if (*det < eps)
        error("Jacobi-Matrix nicht positiv definit (Element %d)", ie);

    xjj = xj[0][0];                              /* Inverse der Jacobi-Matrix */
    xj[0][0] =  xj[1][1] / *det;
    xj[0][1] = -xj[0][1] / *det;
    xj[1][0] = -xj[1][0] / *det;
    xj[1][1] =       xjj / *det;

    for (k = 0; k < 4; ++k) {             /* Verzerrungs-Verschiebungs-Matrix */
        k1 = k * 2;
        k2 = k1 + 1;
        d[0][k1] = xj[0][0] * p[0][k] + xj[0][1] * p[1][k];
        d[0][k2] = 0.0;
        d[1][k1] = 0.0;
        d[1][k2] = xj[1][0] * p[0][k] + xj[1][1] * p[1][k];
        d[2][k1] = d[1][k2];
        d[2][k2] = d[0][k1];
    }

    if (ic == 4) {                          /* rotationssymmetrisches Element */
        h[0] = 0.25 * rp * sp;     /* Ansatzfunktionen fuer r und s auswerten */
        h[1] = 0.25 * rm * sp;
        h[2] = 0.25 * rm * sm;
        h[3] = 0.25 * rp * sm;

        *rad = h[0] * x1[0] +            /* Radius (Y-Achse = Symmetrieachse) */
               h[1] * x2[0] +
               h[2] * x3[0] +
               h[3] * x4[0];
        if (*rad < eps) error("nichtpositiver Radius (Element %d)", ie);

        for (k = 0; k < 4; ++k) {         /* Verzerrungs-Verschiebungs-Matrix */
            d[3][k * 2    ] = h[k] / *rad;
            d[3][k * 2 + 1] = 0.0;
        }
    }
}

/*\
|*|     Elementsteifigkeitsmatrix in die reduzierte Gesamtsteifigkeitsmatrix
|*|     einfuegen
|*|
|*|     Input:
|*|
|*|         m23            = 2 = zweidimensionales Modell
|*|                        = 3 = dreidimensionales Modell 
|*|         nne            = Anzahl der Knoten des Elements
|*|         nele[nne]      = Knotennummern des Elements
|*|         id[nn][m23]    = Freiheitsgradnummern
|*|         s[nd*(nd+1)/2] = Elementsteifigkeitsmatrix
|*|                          (Symmetrie ausgenutzt, nd = m23*nne)
|*|         mda[neq+1]     = Adressen der Diagonalelemente
|*|
|*|     Output:
|*|
|*|         a[nwk]         = reduzierte Gesamtsteifigkeitsmatrix
|*|
|*|         | s[0] s[1]  s[2]      ...  |      | a[0] a[2] a[5] ...  |
|*|         |      s[nd] s[nd+1]   ...  |      |      a[1] a[4] ...  |
|*|     s = |            s[2*nd-1] ...  |, a = |           a[3] ...  |
|*|         |                      ...  |      |                ...  |
\*/
void addele(int m23, int nne, int *nele, int **id, double *s, int *mda,
    double *a)
{
    int i, ii, i1, i2, j, jj, j1, j2, ij, ndi = 0, ks;

    for (i1 = 0; i1 < nne; ++i1) {
        for (i2 = 0; i2 < m23; ++i2) {
            i = i1 * m23 + i2 + 1;
            ii = id[nele[i1] - 1][i2];
            if (ii > 0) {
                ks = i;
                for (j1 = 0; j1 < nne; ++j1) {
                    for (j2 = 0; j2 < m23; ++j2) {
                        j = j1 * m23 + j2 + 1;
                        jj = id[nele[j1] - 1][j2];
                        if (jj > 0 && (ij = ii - jj) >= 0)
                            a[mda[ii - 1] + ij] +=
                                s[((j >= i) ? j + ndi : ks) - 1];
                        ks = ks + m23 * nne - j;
                    }
                }
            }
            ndi = ndi + m23 * nne - i;
        }

    }
}

/*\
|*|     Lineares Gleichungssystem mit symmetrischer Koeffizientenmatrix mit
|*|     Hilfe des Gaussschen Algorithmus loesen
|*|
|*|     Input:
|*|
|*|         ktr        = 0 = Komplette Loesung
|*|                    = 1 = Dreieckszerlegung der Koeffizientenmatrix
|*|                    = 2 = Ermittlung des Loesungsvektors fuer eine rechte
|*|                          Seite
|*|         neq        = Anzahl der Freiheitsgrade
|*|         mda[neq+1] = Adressen der Diagonalelemente
|*|         eps        = Genauigkeit
|*|         a[nwk]     = reduzierte Koeffizientenmatrix bzw. deren
|*|                      Dreieckszerlegung
|*|         v[neq]     = rechte Seite
|*|
|*|     Output:
|*|
|*|         a[nwk]     = Dreieckszerlegung der Koeffizientenmatrix
|*|         v[neq]     = Loesungsvektor
|*|
|*|     Die hier gewaehlte Variation des Gaussschen Verfahrens nach Crout und
|*|     Banachiewicz hat nur zur Voraussetzung, dass keine Komponente der
|*|     Diagonale waehrend der Dreieckszerlegung einen Wert Null oder nahe Null
|*|     erhaelt und waere damit auch fuer nicht positiv definite Matrizen
|*|     anwendbar.
\*/
void colsol(int ktr, int neq, int *mda, const double eps, double *a, double *v)
{
    int n, j, k, l, kn, kl, ku, kh, ic, klt, ki, nd, kk;
    double b, c;

    if (ktr != 2) {              /* Dreieckszerlegung der Koeffizientenmatrix */
        for (n = 0; n < neq; ++n) {
            kn = mda[n];                                   /* Diagonalelement */
            kl = kn + 1;              /* unterstes Element der aktiven Spalte */
            ku = mda[n + 1] - 1;       /* oberstes Element der aktiven Spalte */
            kh = ku - kl;                                 /* Spaltenhoehe - 1 */
            if (kh >= 0) {
                if (kh) {
                    k = n - kh;
                    ic = 0, klt = ku;
                    for (j = 0; j < kh; ++j) {
                        ++ic, --klt;
                        ki = mda[k];
                        nd = mda[k + 1] - ki - 1;
                        if (nd > 0) {
                            kk = (ic < nd) ? ic : nd;
                            c = 0.0;
                            for (l = 1; l <= kk; ++l)
                                c += a[ki + l] * a[klt + l];
                            a[klt] -= c;
                        }
                        ++k;
                    }
                }
                b = 0.0;
                for (kk = kl, k = n; kk <= ku; ++kk) {
                    c = a[kk] / a[mda[--k]];
                    b += c * a[kk];
                    a[kk] = c;
                }
                a[kn] -= b;
            }
            if (a[kn] < eps) error("Koeffizientenmatrix nicht positiv definit,"
                " Pivot[%d] = %g", n, a[kn]);
        }
    }

    if (ktr != 1) {                            /* Reduktion der rechten Seite */
        for (n = 0; n < neq; ++n) {
            kl = mda[n] + 1;
            ku = mda[n + 1] - 1;
            if (ku < kl) continue;
            c = 0.0;
            for (kk = kl, k = n; kk <= ku; ++kk) c += a[kk] * v[--k];
            v[n] -= c;
        }
                                      /* Division durch Diagonalkoeffizienten */
        for (n = 0; n < neq; ++n) v[n] /= a[mda[n]];

        for (l = 1, n = neq; l < neq; ++l) {             /* Ruecksubstitution */
            ku = mda[n] - 1;
            kl = mda[--n] + 1;
            if (ku < kl) continue;
            for (kk = kl, k = n; kk <= ku; ++kk) v[--k] -= a[kk] * v[n];
        }
    }
}

/*\
|*|     Verschiebungen ausgeben
|*|
|*|     Input:
|*|
|*|         m23         = 2 = zweidimensionales Modell
|*|                     = 3 = dreidimensionales Modell
|*|         nn          = Anzahl der Knoten
|*|         neq         = Anzahl der Freiheitsgrade
|*|         id[nn][m23] = Freiheitsgradnummern 
|*|         v[neq]      = Verschiebungsvektor
|*|         out         = Ausgabestrom
\*/
void displacements(int m23, int nn, int neq, int **id, double *v, FILE *out)
{
    int i, j, k;

    print(out, "\nKnoten Verschiebungen\n");
    for (i = 0; i < nn; ++i) {
        print(out, "%d  ", i + 1);
        for (j = 0; j < m23; ++j) {
            k = id[i][j] - 1;
            print(out, " %g", (k >= 0 && k < neq) ? v[k] : 0.0);
        }
        print(out, "\n");
    }
}