输入是结点力，输出是位移

snow.cpp

#include "include.hpp"
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
using std::cout;
using std::endl;
using std::string;
using std::ifstream;
using std::ofstream;
using std::vector;
bool read(size_t& nx1, size_t& ny1, size_t& nx2, size_t& ny2, istream& in, size_t jump);
bool read(size_t& nx1, size_t& ny1, size_t& nx2, size_t& ny2, double& q, istream& in);
bool read(double& e, double& mu, istream& in);
int main()
{
    size_t nx1, ny1, nx2, ny2;
    Mesh m;
    ifstream in;

    //读取结点力
    in.open("input.txt");
    read(nx1, ny1, nx2, ny2, in, 0);
    in.close();

    //解方程，Jacobi行列式
    double dy1 = 10.0 / ny1;
    size_t i, j;
    for (i = 0; i != ny1; ++i) {
        double y0 = i * dy1;
        double y1 = y0 + dy1;
        double dx0 = (14.0 - y0 / 2.5) / nx1;
        double dx1 = (14.0 - y1 / 2.5) / nx1;
        for (j = 0; j != nx1; ++j) {
            double x00 = y0 / 2.5 + j * dx0;
            double x10 = x00 + dx0;
            double x01 = y1 / 2.5 + j * dx1;
            double x11 = x01 + dx1;
            m.addnode(i * (nx1 + 1) + j, Point(x00, y0));
            m.addnode(i * (nx1 + 1) + j + 1, Point(x10, y0));
            m.addnode((i + 1) * (nx1 + 1) + j + 1, Point(x11, y1));
            m.addnode((i + 1) * (nx1 + 1) + j, Point(x01, y1));
            m.addunit(QuadS(i * (nx1 + 1) + j, i * (nx1 + 1) + j + 1, (i + 1) * (nx1 + 1) + j + 1, (i + 1) * (nx1 + 1) + j), 0);
        }
    }
    NodeData<double> rx, ry;
    NodeData<bool> bx, by;
    Mesh::node_i ni;
    for (ni = m.nodebegin(); ni != m.nodeend(); ++ni) {
        bx.add(0);
        rx.add(0);
        by.add(0);
        ry.add(0);
    }
    for (j = 0; j != nx1 + 1; ++j) {
        bx.set(j, 1);
        by.set(j, 1);
    }
    int a, b;
    double q;
    
    //读取结点力
    in.open("input.txt");
    read(nx1, ny1, nx2, ny2, q, in);
    in.close();

    size_t end = m.nodenum() - 1;
    NodeData<double> loadx;
    double f = 540 / (ny1 ^ 2) + 540 * (2 * ny1 - 2 * j - 1) / (ny1 ^ 2);
    for (ni = m.nodebegin(); ni != m.nodeend(); ++ni) {
        loadx.add(ni->first, 0);
    }
    loadx.set(0, 0.928 * (540 * ny1 - 180) / (ny1 ^ 2));
    loadx.set(end - nx1, 0.928 * 180 / ((ny1) ^ 2));
    for (j = 1; j != ny1 - 1; j++) {
        a = (nx1 + 1) * j;
        loadx.set(a, 0.928 * f);
    }

    NodeData<double> loady;
    for (ni = m.nodebegin(); ni != m.nodeend(); ++ni) {
        loady.add(ni->first, 0);
    }
    loady.set(0, 0.37 * (540 * ny1 - 180) / (ny1 ^ 2));
    loady.set(end - nx1, 0.37 * 180 / ((ny1) ^ 2));
    for (j = 1; j != ny1 - 1; j++) {
        b = (nx1 + 1) * j;
        loady.set(b, -0.37 * f);
    }
    for (ni = m.nodebegin(); ni != m.nodeend(); ++ni) {
        loady.add(ni->first, 0);
    }

    //读取结点力
   	double e, mu;
    in.open("input.txt");
    read(e, mu, in);
    in.close();
    Material mat;
    vector<double> datau;
    datau.push_back(e);
    datau.push_back(mu);
    mat.data[0] = datau;
    cout << "组装总刚" << endl;
    Tiff K;
    cout << "预组装";
    K.premake(m);
    cout << "\t完成" << endl;
    cout << "********************************" << endl;
    cout << "计算总刚";
    K.make(m, mat);
    cout << "\t完成" << endl;
    cout << "********************************" << endl;
    cout << "引入载荷、约束";
    K.posmake(bx, rx, loadx, by, ry, loady);
    bx.clear();
    by.clear();
    rx.clear();
    ry.clear();
    cout << "\t完成" << endl;


    //开始解方程
    cout << "********************************" << endl
         << endl;
    cout << "求解方程";
    K.solve(loadx, loady);
    K.claer();
    cout << "\t完成" << endl;
    cout << "********************************" << endl
         << endl;
    cout << "计算应力";
    vector<NodeData<double> > s = gets(loadx, loady, m, mat);
    cout << "\t完成" << endl;

    //输出结果到文件
    cout << "********************************" << endl
         << endl;
    cout << "输出文件";
    vector<string> datahead;
    datahead.push_back("U");
    datahead.push_back("V");
    datahead.push_back("Fx");
    datahead.push_back("Fy");
    datahead.push_back("Fxy");
    vector<NodeData<double> > output;
    output.push_back(loadx);
    output.push_back(loady);
    output.insert(output.end(), s.begin(), s.end());
    ofstream out;
    out.open("output.dat");
    m.show_head("Output File", datahead, out);
    m.show_data(output, out);
    m.show_unit(out);
    out.close();
    cout << "\t完成" << endl;
    cout << "********************************" << endl
         << endl;
    return 0;
}


/*********************************
* 读取字节
*********************************/
bool read(size_t& nx1, size_t& ny1, size_t& nx2, size_t& ny2, istream& in, size_t jump)
{
    char temp;
    for (size_t i = 0; i != jump; ++i) {
        do {
            in.get(temp);
        } while (temp != '=');
    }
    do {
        in.get(temp);
    } while (temp != '=');
    in >> nx1;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> ny1;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> nx2;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> ny2;
    return true;
}


/*********************************
* 读取数字
*********************************/
bool read(double& e, double& mu, istream& in)
{
    char temp;
    size_t i;
    for (i = 0; i != 5; ++i) {
        do {
            in.get(temp);
        } while (temp != '=');
    }
    in >> e;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> mu;
    return true;
}

/*********************************
* 读取字节
*********************************/
bool read(size_t& nx1, size_t& ny1, size_t& nx2, size_t& ny2, double& q, istream& in)
{
    char temp;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> nx1;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> ny1;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> nx2;
    do {
        in.get(temp);
    } while (temp != '=');
    in >> ny2;
    for (size_t i = 0; i != 3; ++i) {
        do {
            in.get(temp);
        } while (temp != '=');
    }
    in >> q;
    return true;
}