fem

//#pragma warning(disable:4996)
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <locale.h>

#include <iostream>
#include <fstream>

#include <algorithm>
#include <vector>
#include <list>
#include <set>
#include <tuple>
#include <utility> // pair

using namespace std;

const double BIG = 1e8;

const int Gx[8][8] = {
    { 4, -4, 2, -2, 2, -2, 1, -1 },
    { -4, 4, -2, 2, -2, 2, -1, 1 },
    { 2, -2, 4, -4, 1, -1, 2, -2 },
    { -2, 2, -4, 4, -1, 1, -2, 2 },
    { 2, -2, 1, -1, 4, -4, 2, -2 },
    { -2, 2, -1, 1, -4, 4, -2, 2 },
    { 1, -1, 2, -2, 2, -2, 4, -4 },
    { -1, 1, -2, 2, -2, 2, -4, 4 },
};

const int Gy[8][8] = {
    { 4, 2, -4, -2, 2, 1, -2, -1 },
    { 2, 4, -2, -4, 1, 2, -1, -2 },
    { -4, -2, 4, 2, -2, -1, 2, 1 },
    { -2, -4, 2, 4, -1, -2, 1, 2 },
    { 2, 1, -2, -1, 4, 2, -4, -2 },
    { 1, 2, -1, -2, 2, 4, -2, -4 },
    { -2, -1, 2, 1, -4, -2, 4, 2 },
    { -1, -2, 1, 2, -2, -4, 2, 4 },
};

const int Gz[8][8] = {
    { 4, 2, 2, 1, -4, -2, -2, -1 },
    { 2, 4, 1, 2, -2, -4, -1, -2 },
    { 2, 1, 4, 2, -2, -1, -4, -2 },
    { 1, 2, 2, 4, -1, -2, -2, -4 },
    { -4, -2, -2, -1, 4, 2, 2, 1 },
    { -2, -4, -1, -2, 2, 4, 1, 2 },
    { -2, -1, -4, -2, 2, 1, 4, 2 },
    { -1, -2, -2, -4, 1, 2, 2, 4 },
};

const int M[8][8] = {
    { 8, 4, 4, 2, 4, 2, 2, 1 },
    { 4, 8, 2, 4, 2, 4, 1, 2 },
    { 4, 2, 8, 4, 2, 1, 4, 2 },
    { 2, 4, 4, 8, 1, 2, 2, 4 },
    { 4, 2, 2, 1, 8, 4, 4, 2 },
    { 2, 4, 1, 2, 4, 8, 2, 4 },
    { 2, 1, 4, 2, 4, 2, 8, 4 },
    { 1, 2, 2, 4, 2, 4, 4, 8 },
};

const int C[4][4] = {
    { 4, 2, 2, 1 },
    { 2, 4, 1, 2 },
    { 2, 1, 4, 2 },
    { 1, 2, 2, 4 },
};

struct node
{
    double x, y, z;
};

struct u_in_node
{
    int number;
    double value;
};

struct area
{
    vector<double> values;
    vector<int> inodes;
};

struct element
{
    double hx, hy, hz, lambda;
    vector<int> inodes;
};

double beta;
double gamma;

double get_h(vector<double> values)
{
    for (int i = 0; i < values.size(); i++)
        for (int j = i + 1; j < values.size(); j++)
        {
            double diff = values[i] - values[j];
            if (diff != 0)
                return abs(diff);
        }
}

void get_all_h(vector<node> nodes, vector<int> inodes, double& hx, double& hy, double& hz)
{
    vector<double> values;
    for (int inode : inodes)
        values.push_back(nodes[inode].x);
    hx = get_h(values);

    values.clear();
    for (int inode : inodes)
        values.push_back(nodes[inode].y);
    hy = get_h(values);

    values.clear();
    for (int inode : inodes)
        values.push_back(nodes[inode].z);
    hz = get_h(values);
}

int input(
    vector<node>& nodes,
    vector<element>& elements,
    vector<double>& f,
    vector<u_in_node>& conditions1,
    vector<area>& conditions2,
    vector<area>& conditions3
)
{
    ifstream file("nodes.txt");
    int n;
    file >> n;
    file >> beta;
    file >> gamma;
    double x, y, z;
    for (int i = 0; i < n; i++)
    {
        file >> x >> y >> z;
        node _node;
        _node.x = x;
        _node.y = y;
        _node.z = z;
        nodes.push_back(_node);
    }
    file.close();

    ifstream file2("elements.txt");
    int elem_count;
    file2 >> elem_count;
    for (int i = 0; i < elem_count; i++)
    {
        element new_elem;
        int temp;
        for (int j = 0; j < 8; j++)
        {
            file2 >> temp;
            new_elem.inodes.push_back(temp - 1);
        }
        file2 >> temp;
        new_elem.lambda = temp;
        sort(new_elem.inodes.begin(), new_elem.inodes.end());
        get_all_h(nodes, new_elem.inodes, new_elem.hx, new_elem.hy, new_elem.hz);
        elements.push_back(new_elem);
    }
    file2.close();

    ifstream file3("f.txt");
    for (int i = 0; i < n; i++)
    {
        file3 >> x;
        f.push_back(x);
    }
    file3.close();

    ifstream file4("cond1.txt");
    int nn;
    file4 >> nn;
    for (int i = 0; i < nn; i++)
    {
        u_in_node _u;
        int temp;
        file4 >> temp;
        _u.number = temp - 1;
        file4 >> _u.value;
        conditions1.push_back(_u);
    }
    file4.close();

    ifstream file5("cond2.txt");
    file5 >> nn;
    double temp2;
    for (int i = 0; i < nn; i++)
    {
        area _area;
        int temp;
        for (int j = 0; j < 4; j++)
        {
            file5 >> temp;
            _area.inodes.push_back(temp - 1);
        }
        double temp2;
        for (int j = 0; j < 4; j++)
        {
            file5 >> temp2;
            _area.values.push_back(temp2);
        }
        conditions2.push_back(_area);
    }
    file5.close();

    ifstream file6("cond3.txt");
    file6 >> nn;
    for (int i = 0; i < nn; i++)
    {
        area _area;
        _area.values.push_back(beta);
        int temp;
        for (int j = 0; j < 4; j++)
        {
            file6 >> temp;
            _area.inodes.push_back(temp - 1);
        }
        double temp2;
        for (int j = 0; j < 4; j++)
        {
            file6 >> temp2;
            _area.values.push_back(temp2);
        }
        conditions3.push_back(_area);
    }
    file6.close();

    return n;
}

void local_matrix(int n, double lambda, double hx, double hy, double hz, vector<double> f, double A[8][8], double b[8])
{
    for (int i = 0; i < 8; i++)
        for (int j = 0; j < 8; j++)
        {
            A[i][j] = lambda * hx * hy * hz / 36 * (Gx[i][j] / (hx * hx) + Gy[i][j] / (hy * hy) + Gz[i][j] / (hz * hz));
            A[i][j] += gamma * hx * hy * hz / 216 * M[i][j];
        }

    for (int i = 0; i < 8; i++)
    {
        double sum = 0;
        for (int j = 0; j < 8; j++)
            sum += M[i][j] * f[j];
        b[i] = gamma * hx * hy * hz * sum / 216;
    }
}

void portrait(vector<element> elements, int n, vector<int>& ia, vector<int>& ja)
{
    vector<set<int>> _set;
    _set.resize(n);

    for (int i = 0; i < elements.size(); i++)
        for (int el : elements[i].inodes)
            for (int node : elements[i].inodes)
                if (el > node)
                    _set[el].insert(node);

    for (int i = 0; i < _set.size(); i++)
    {
        printf("%d : ", i + 1);
        for (int a : _set[i])
            printf("%d ", a + 1);
        cout << endl;
    }

    ia.push_back(0);

    for (int i = 1; i < _set.size(); i++)
    {
        ia.push_back(ia.back() + _set[i].size());
        for (int node : _set[i])
            ja.push_back(node);
    }

    for (int i = 0; i < ia.size(); i++)
    {
        if (ia[i] == 0)
            ia[i] = 1;
        else break;
    }

    cout << "\nportrait start\n";

    for (int i = 0; i < ia.size(); i++)
        printf("%d ", ia[i]);

    cout << endl;

    for (int i = 0; i < ja.size(); i++)
        printf("%d ", ja[i] + 1);

    cout << "\nportrait end\n\n";
}

void matrix_build(int n, vector<node> nodes, vector<element> elements, vector<double> f, double** A, double* b)
{
    /*double** A = new double* [n];
    for (int i = 0; i < n; i++)
        A[i] = new double[n];
    double* b = new double[n];*/

    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            A[i][j] = 0;
        b[i] = 0;
    }

    double _A[8][8];
    double _b[8];

    for (int i = 0; i < elements.size(); i++)
    {
        local_matrix(n, elements[i].lambda, elements[i].hx, elements[i].hy, elements[i].hz, f, _A, _b);
        for (int j = 0; j < 8; j++)
        {

            for (int k = 0; k < 8; k++)
            {
                A[elements[i].inodes[j]][elements[i].inodes[k]] += _A[j][k];
                //printf("%d -> %d\n", elements[i].inodes[j], j);
                //printf("%d -> %d\n", elements[i].inodes[k], k);
            }
            b[elements[i].inodes[j]] += _b[j];
        }

        printf("\n\n");
        for (int i = 0; i < 8; i++)
        {
            for (int j = 0; j < 8; j++)
                printf("%5.2f ", _A[i][j]);
            cout << endl;
        }
        printf("\n\n");
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
                printf("%5.2f ", A[i][j]);
            cout << endl;
        }
    }
}

void conditions(vector<node> nodes, vector<u_in_node> conditions1, vector<area> conditions2, vector<area> conditions3, double** A, double* b)
{
    cout << "\n\nconditions\n\n";

    for (area _a : conditions3)
    {
        double hx = abs(nodes[_a.inodes[0]].x - nodes[_a.inodes[3]].x);
        double hy = abs(nodes[_a.inodes[0]].y - nodes[_a.inodes[3]].y);
        double hz = abs(nodes[_a.inodes[0]].z - nodes[_a.inodes[3]].z);
        hx = (hx == 0) ? 1 : hx;
        hy = (hy == 0) ? 1 : hy;
        hz = (hz == 0) ? 1 : hz;

        double _C[4];
        for (int i = 0; i < 4; i++)
        {
            _C[i] = 0;
            for (int j = 0; j < 4; j++)
                _C[i] += C[i][j] * _a.values[j];
        }

        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
                A[_a.inodes[i]][_a.inodes[j]] += beta * hx * hy * hz / 36 * C[i][j];
            b[_a.inodes[i]] += beta * hx * hy * hz / 36 * _C[i];
        }

        printf("\n\n");
        for (int i = 0; i < nodes.size(); i++)
        {
            for (int j = 0; j < nodes.size(); j++)
                printf("%5.2f ", A[i][j]);
            cout << endl;
        }
        printf("\n\n");
        for (int i = 0; i < nodes.size(); i++)
            printf("%5.2f\n", b[i]);
    }

    for (area _a : conditions2)
    {
        double hx = abs(nodes[_a.inodes[0]].x - nodes[_a.inodes[3]].x);
        double hy = abs(nodes[_a.inodes[0]].y - nodes[_a.inodes[3]].y);
        double hz = abs(nodes[_a.inodes[0]].z - nodes[_a.inodes[3]].z);
        hx = (hx == 0) ? 1 : hx;
        hy = (hy == 0) ? 1 : hy;
        hz = (hz == 0) ? 1 : hz;

        double _C[4];
        for (int i = 0; i < 4; i++)
            for (int j = 0; j < 4; j++)
                b[_a.inodes[i]] += beta * hx * hy * hz / 36 * C[i][j] * _a.values[j];

        printf("\n\n");
        for (int i = 0; i < nodes.size(); i++)
            printf("%5.2f\n", b[i]);
    }

    for (u_in_node _u : conditions1)
    {
        A[_u.number][_u.number] = BIG;
        b[_u.number] = BIG * _u.value;
    }
}

int main()
{
    vector<node> nodes;
    vector<element> elements;
    vector<double> f;

    vector<u_in_node> conditions1;
    vector<area> conditions2;
    vector<area> conditions3;

    vector<int> ia = { 0 };
    vector<int> ja;

    /*double _A[8][8];
    double _b[8];*/

    int n = input(nodes, elements, f, conditions1, conditions2, conditions3);

    //portrait(elements, n, ia, ja);
    //matrix_build(n, nodes, elements, f, ia, ja);

    double** A = new double* [n];
    for (int i = 0; i < n; i++)
        A[i] = new double[n];
    double* b = new double[n];

    matrix_build(n, nodes, elements, f, A, b);

    conditions(nodes, conditions1, conditions2, conditions3, A, b);

    return 0;
}