initialize

.gitignore

@@ -0,0 +1,8 @@
+# vscode
+.vscode
+
+# compile
+build
+
+# output
+output

cpp/CMakeLists.txt

@@ -0,0 +1,26 @@
+cmake_minimum_required(VERSION 3.9)
+project(fem3d)
+
+# project source files
+file(GLOB SRCFILES
+    "src/*.cpp"
+    "src/LinSysSolver/*.cpp"
+)
+
+FOREACH(item ${SRCFILES})
+    IF(${item} MATCHES "main.cpp")
+        LIST(REMOVE_ITEM SRCFILES ${item})
+    ENDIF(${item} MATCHES "main.cpp")
+ENDFOREACH(item)
+
+add_library(${PROJECT_NAME}_dev ${SRCFILES})
+
+target_include_directories(${PROJECT_NAME}_dev PUBLIC
+    "src"
+    "src/LinSysSolver"
+)
+
+add_executable(${PROJECT_NAME}_bin "src/main.cpp")
+target_link_libraries(${PROJECT_NAME}_bin PUBLIC ${PROJECT_NAME}_dev)
+# add path of Eigen header file to include directories
+target_include_directories(${PROJECT_NAME}_dev PUBLIC "/usr/local/include/eigen-3.4.0")

cpp/Readme.md

@@ -0,0 +1,35 @@
+# Readme
+
+### 1. 依赖
+
+* `Eigen`
+
+### 2. 运行
+
+```shell
+mkdir build
+cd build
+cmake ..
+make
+```
+
+* `Eigen`: 运行前，需要先在`CMakeLists.txt`的第26行修改`Eigen`的目录路径
+* 模型的杨氏模量、泊松比、热膨胀系数、参考温度`Tref`：通过`src/Mesh.cpp`的17~20行分别设定
+* 固定点`DBC`：通过修改`src/Mesh.cpp`的`setDBC`函数设定
+
+### 3. 输入
+
+* `input/cube.txt`：自定义格式，四面体网格信息；
+
+  ​								  第一、二行分别是顶点数量`n`与四面体单元数量`m`，
+
+  ​								  从第三行开始是`n`行顶点坐标与`m`行四面体顶点下标，下标从1开始
+
+* `input/cube_T.txt`：热传导仿真结果，`n`行温度数值，单位开尔文（K）
+
+### 4. 输出
+
+* `output/U.txt`：四面体网格每个顶点的位移，
+
+  ​			   		  	`n*3`行，第$3i-2$，$3i-1$，$3i$行表示第$i$个节点在$x$，$y$，$z$轴上的位移，$1 \leq i \leq n$
+

cpp/src/LinSysSolver/EigenLibSolver.cpp

@@ -0,0 +1,77 @@
+//
+// EigenLibSolver.cpp
+//
+// Created by Wei Chen on 9/2/21
+//
+
+#include "EigenLibSolver.hpp"
+
+namespace fem3d{
+
+void EigenLibSolver::set_pattern(const std::vector<std::set<int>>& vNeighbor){
+    Base::set_pattern(vNeighbor);
+
+    coefMtr.resize(Base::numRows, Base::numRows);
+    coefMtr.reserve(Base::nnz);
+
+    memcpy(coefMtr.innerIndexPtr(), Base::ja.data(), Base::ja.size()*sizeof(Base::ja[0]));
+    memcpy(coefMtr.outerIndexPtr(), Base::ia.data(), Base::ia.size()*sizeof(Base::ia[0]));
+}
+
+void EigenLibSolver::analyze_pattern(void){
+    simplicialLDLT.analyzePattern(coefMtr);
+    assert(simplicialLDLT.info() == Eigen::Success);
+}
+
+void EigenLibSolver::factorize(void){
+    simplicialLDLT.factorize(coefMtr);
+    assert(simplicialLDLT.info() == Eigen::Success);
+}
+
+void EigenLibSolver::solve(const Eigen::VectorXd& rhs, Eigen::VectorXd& res){
+    res = simplicialLDLT.solve(rhs);
+    assert(simplicialLDLT.info() == Eigen::Success);
+}
+
+void EigenLibSolver::setCoeff(int rowI, int colI, double val){
+    assert(rowI < Base::numRows);
+    const auto finder = Base::IJ2aI[rowI].find(colI);
+    assert(finder != Base::IJ2aI[rowI].end());
+    Base::a[finder->second] = val;
+    coefMtr.valuePtr()[finder->second] = val;
+}
+
+void EigenLibSolver::addCoeff(int rowI, int colI, double val){
+    assert(rowI < Base::numRows);
+    const auto finder = Base::IJ2aI[rowI].find(colI);
+    assert(finder != Base::IJ2aI[rowI].end());
+    Base::a[finder->second] += val;
+    coefMtr.valuePtr()[finder->second] += val;
+}
+
+void EigenLibSolver::setZero(void){
+    Base::setZero();
+    memcpy(coefMtr.valuePtr(), Base::a.data(), Base::a.size()*sizeof(Base::a[0]));
+}
+
+void EigenLibSolver::setUnitRow(int rowI){
+    assert(rowI < Base::numRows);
+    for(const auto& colI : Base::IJ2aI[rowI]){
+        double tmp = (rowI == colI.first);
+        Base::a[colI.second] = tmp;
+        coefMtr.valuePtr()[colI.second] = tmp;
+    }
+}
+
+void EigenLibSolver::setZeroCol(int colI){
+    assert(colI < Base::numRows);
+    for(int rowI=0; rowI<Base::numRows; ++rowI){
+        const auto finder = Base::IJ2aI[rowI].find(colI);
+        if(finder != Base::IJ2aI[rowI].end()){
+            Base::a[finder->second] = 0.0;
+            coefMtr.valuePtr()[finder->second] = 0.0;
+        }
+    }
+}
+
+} // namespace

cpp/src/LinSysSolver/EigenLibSolver.hpp

@@ -0,0 +1,37 @@
+//
+// EigenLibSolver.hpp
+//
+// Created by Wei Chen on 9/2/21
+//
+
+#ifndef EigenLibSolver_hpp
+#define EigenLibSolver_hpp
+
+#include "LinSysSolver.hpp"
+
+namespace fem3d{
+
+class EigenLibSolver : public LinSysSolver{
+    typedef LinSysSolver Base;
+
+protected:
+    Eigen::SparseMatrix<double> coefMtr;
+    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> simplicialLDLT;
+
+public:
+    virtual void set_pattern(const std::vector<std::set<int>>& vNeighbor);
+
+    virtual void analyze_pattern(void);
+    virtual void factorize(void);
+    virtual void solve(const Eigen::VectorXd& rhs, Eigen::VectorXd& res);
+
+    virtual void setCoeff(int rowI, int colI, double val);
+    virtual void addCoeff(int rowI, int colI, double val);
+    virtual void setZero(void);
+    virtual void setUnitRow(int rowI);
+    virtual void setZeroCol(int colI);
+};
+
+} // namespace
+
+#endif // EigenLibSolver_hpp

cpp/src/LinSysSolver/LinSysSolver.hpp

@@ -0,0 +1,113 @@
+//
+// LinSysSolver.hpp
+//
+// Created by Wei Chen on 9/2/21
+//
+
+#ifndef LinSysSolver_hpp
+#define LinSysSolver_hpp
+
+#include <Eigen/Eigen>
+#include <iostream>
+#include <set>
+#include <map>
+#include <vector>
+
+namespace fem3d{
+
+class LinSysSolver{
+protected:
+    int numRows, nnz;
+    Eigen::VectorXi ia, ja;
+    std::vector<std::map<int, int>> IJ2aI;
+    Eigen::VectorXd a;
+
+public:
+    virtual ~LinSysSolver(void){};
+
+    virtual void set_pattern(const std::vector<std::set<int>>& vNeighbor){
+        int nvet = static_cast<int>(vNeighbor.size());
+        numRows = nvet * 3;
+        nnz = 0;
+        for(int vI=0; vI<nvet; ++vI){
+            nnz += static_cast<int>(vNeighbor[vI].size()) * 9;
+        }
+
+        ia.setZero(numRows+1);
+        ja.setZero(nnz);
+        IJ2aI.resize(0);
+        IJ2aI.resize(numRows);
+        a.setZero(nnz);
+
+        for(int vI=0; vI<nvet; ++vI){
+            int num = static_cast<int>(vNeighbor[vI].size()) * 3;
+            for(int dof=0; dof<3; ++dof){
+                int row = vI * 3 + dof;
+                ia(row+1) = ia(row) + num;
+                int cnt = 0;
+                for(const auto& nbVI : vNeighbor[vI]){
+                    ja(ia(row)+cnt) = nbVI * 3;
+                    ja(ia(row)+cnt+1) = nbVI * 3 + 1;
+                    ja(ia(row)+cnt+2) = nbVI * 3 + 2;
+                    IJ2aI[row][nbVI*3] = ia(row) + cnt;
+                    IJ2aI[row][nbVI*3+1] = ia(row) + cnt + 1;
+                    IJ2aI[row][nbVI*3+2] = ia(row) + cnt + 2;
+                    cnt += 3;
+                }
+            }
+        }
+    }
+
+    virtual void analyze_pattern(void) = 0;
+    virtual void factorize(void) = 0;
+    virtual void solve(const Eigen::VectorXd& rhs, Eigen::VectorXd& res) = 0;
+
+    virtual void setCoeff(int rowI, int colI, double val){
+        assert(rowI<numRows);
+        const auto finder = IJ2aI[rowI].find(colI);
+        assert(finder!=IJ2aI[rowI].end());
+        a[finder->second] = val;
+    }
+
+    virtual void addCoeff(int rowI, int colI, double val){
+        assert(rowI<numRows);
+        const auto finder = IJ2aI[rowI].find(colI);
+        assert(finder!=IJ2aI[rowI].end());
+        a[finder->second] += val;
+    }
+
+    virtual void setZero(void){
+        a.setZero();
+    }
+
+    virtual void multiply(const Eigen::VectorXd& x, Eigen::VectorXd& Ax){
+        assert(x.size()==numRows);
+        Ax.setZero(numRows);
+        for(int rowI=0; rowI<numRows; ++rowI){
+            for(const auto& colI : IJ2aI[rowI]){
+                Ax(rowI) += a[colI.second] * x(colI.first);
+            }
+        }
+    }
+
+    virtual void setUnitRow(int rowI){
+        assert(rowI<numRows);
+        for(const auto& colI : IJ2aI[rowI]){
+            a[colI.second] = (rowI==colI.first);
+        }
+    }
+
+    virtual void setZeroCol(int colI){
+        assert(colI<numRows);
+        for(int rowI=0; rowI<numRows; ++rowI){
+            const auto finder = IJ2aI[rowI].find(colI);
+            if(finder != IJ2aI[rowI].end()){
+                a[finder->second] = 0.0;
+            }
+        }
+    }
+};
+
+} // namespace
+
+#endif // LinSysSolver_hpp

cpp/src/Mesh.cpp

@@ -0,0 +1,223 @@
+//
+// Mesh.cpp
+//
+// Created by Wei Chen on 9/1/21
+//
+
+#include "Mesh.hpp"
+#include "Utils.hpp"
+
+#include <iostream>
+#include <cmath>
+#include <cstdlib>
+
+namespace fem3d{
+
+Mesh::Mesh(void){
+    E = 73.0;
+    nu = 0.17;
+    alpha = 5.5e-7;
+    Tref = 293.15;
+    D.resize(6, 6);
+    D << 1.0-nu,     nu,     nu,              0.0,              0.0,              0.0,
+             nu, 1.0-nu,     nu,              0.0,              0.0,              0.0,
+             nu,     nu, 1.0-nu,              0.0,              0.0,              0.0,
+            0.0,    0.0,    0.0, (1.0-2.0*nu)/2.0,              0.0,              0.0,
+            0.0,    0.0,    0.0,              0.0, (1.0-2.0*nu)/2.0,              0.0,
+            0.0,    0.0,    0.0,              0.0,              0.0, (1.0-2.0*nu)/2.0;
+    double tmp = E / (1.0 + nu) / (1.0 - 2.0 * nu);
+    D = D * tmp;
+
+    if(!Utils::readMesh("../input/cube.txt", nvet, nele, vet, ele)){
+        std::cout << "Logging: [fem3d] error on reading vertices" << std::endl;
+        exit(-1);
+    }
+    T.resize(nvet);
+    if(!Utils::readTemp("../input/cube_T.txt", T)){
+        std::cout << "Logging: [fem3d] error on reading temperature" << std::endl;
+        exit(-1);
+    }
+    T = alpha * (T - Eigen::VectorXd::Constant(nvet, Tref));
+    ndof = nvet * 3;
+
+    F.setZero(ndof);
+    U.setZero(ndof);
+
+    computeFeatures();
+    linSysSolver = new EigenLibSolver();
+    linSysSolver->set_pattern(vNeighbor);
+    linSysSolver->analyze_pattern();
+    std::cout << "Logging: [fem3d] mesh constructed" << std::endl;
+}
+
+Mesh::~Mesh(void){
+    delete linSysSolver;
+}
+
+void Mesh::computeFeatures(void){ // compute vNeighbor
+    vNeighbor.resize(0);
+    vNeighbor.resize(nvet);
+
+    for(int vI=0; vI<nvet; ++vI){
+        vNeighbor[vI].insert(vI);
+    }
+
+    for(int eleI=0; eleI<nele; ++eleI){
+        const Eigen::Matrix<int, 1, 4>& eleVInd = ele.row(eleI);
+        for(int i=0; i<4; ++i){
+            for(int j=i+1; j<4; ++j){
+                vNeighbor[eleVInd(i)].insert(eleVInd(j));
+                vNeighbor[eleVInd(j)].insert(eleVInd(i));
+            }
+        }
+    }
+}
+
+void Mesh::computeK(void){ // assembly stiffness matrix K and load F
+    linSysSolver->setZero();
+    
+    for(int eleI=0; eleI<nele; ++eleI){
+        const Eigen::Matrix<int, 1, 4>& eleVInd = ele.row(eleI);
+        Eigen::Matrix<double, 4, 3> Xe;
+        Xe.row(0) = vet.row(eleVInd(0));
+        Xe.row(1) = vet.row(eleVInd(1));
+        Xe.row(2) = vet.row(eleVInd(2));
+        Xe.row(3) = vet.row(eleVInd(3));
+        double Te = (T(eleVInd(0)) + T(eleVInd(1)) + T(eleVInd(2)) +
+                     T(eleVInd(3))) / 4.0;
+
+        int num = 12;
+        Eigen::MatrixXd Ke(num, num);
+        Eigen::VectorXd Fe(num);
+        computeKe(Xe, Te, Ke, Fe);
+
+        Eigen::VectorXi edof(num);
+        edof << eleVInd(0)*3, eleVInd(0)*3+1, eleVInd(0)*3+2,
+                eleVInd(1)*3, eleVInd(1)*3+1, eleVInd(1)*3+2,
+                eleVInd(2)*3, eleVInd(2)*3+1, eleVInd(2)*3+2,
+                eleVInd(3)*3, eleVInd(3)*3+1, eleVInd(3)*3+2;
+                
+        for(int i=0; i<num; ++i){ // assembly K
+            for(int j=0; j<num; ++j){
+                linSysSolver->addCoeff(edof(i), edof(j), Ke(i, j));
+            }
+        }
+
+        for(int i=0; i<num; ++i){ // assembly F
+            F(edof(i)) += Fe(i);
+        }
+    }
+}
+
+// compute element stiffness matrix
+void Mesh::computeKe(const Eigen::Matrix<double, 4, 3>& X, double T,
+                     Eigen::MatrixXd& Ke, Eigen::VectorXd& Fe){
+    Eigen::MatrixXd H(4, 4);
+    H << 1.0, X(0, 0), X(0, 1), X(0, 2),
+         1.0, X(1, 0), X(1, 1), X(1, 2),
+         1.0, X(2, 0), X(2, 1), X(2, 2),
+         1.0, X(3, 0), X(3, 1), X(3, 2);
+    double V6 = H.determinant();
+
+    Eigen::VectorXi rowInd(3);
+    Eigen::VectorXi colInd(3);
+    Eigen::MatrixXd sub(3, 3);
+
+    colInd << 0, 2, 3;
+    rowInd << 1, 2, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double b1 = -sub.determinant();
+    rowInd << 0, 2, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double b2 = sub.determinant();
+    rowInd << 0, 1, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double b3 = -sub.determinant();
+    rowInd << 0, 1, 2;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double b4 = sub.determinant();
+
+    colInd << 0, 1, 3;
+    rowInd << 1, 2, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double c1 = sub.determinant();
+    rowInd << 0, 2, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double c2 = -sub.determinant();
+    rowInd << 0, 1, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double c3 = sub.determinant();
+    rowInd << 0, 1, 2;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double c4 = -sub.determinant();
+
+
+    colInd << 0, 1, 2;
+    rowInd << 1, 2, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double d1 = -sub.determinant();
+    rowInd << 0, 2, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double d2 = sub.determinant();
+    rowInd << 0, 1, 3;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double d3 = -sub.determinant();
+    rowInd << 0, 1, 2;
+    Utils::subMatrix(H, rowInd, colInd, sub);
+    double d4 = sub.determinant();
+
+    Eigen::MatrixXd B(6, 12);
+    B << b1,  0,  0, b2,  0,  0, b3,  0,  0, b4,  0,  0,
+          0, c1,  0,  0, c2,  0,  0, c3,  0,  0, c4,  0,
+          0,  0, d1,  0,  0, d2,  0,  0, d3,  0,  0, d4,
+         c1, b1,  0, c2, b2,  0, c3, b3,  0, c4, b4,  0,
+          0, d1, c1,  0, d2, c2,  0, d3, c3,  0, d4, c4,
+         d1,  0, b1, d2,  0, b2, d3,  0, b3, d4,  0, b4;
+    B = B / V6;
+
+    Ke = V6 / 6.0 * (B.transpose() * D * B);
+
+    Eigen::VectorXd strain(6);
+    strain << T, T, T, 0.0, 0.0, 0.0;
+    Fe = V6 / 6.0 * (B.transpose() * D * strain);
+}
+
+void Mesh::boundaryK(void){
+    for(int i=0; i<DBCV.size(); ++i){
+        int DBCI = DBCV(i);
+        linSysSolver->setZeroCol(DBCI);
+        linSysSolver->setUnitRow(DBCI);
+    }
+}
+
+void Mesh::solve(void){
+    linSysSolver->factorize();
+    linSysSolver->solve(F, U);
+
+    for(int i=0; i<DBCV.size(); ++i){ // boundary DBC dofs
+        U(DBCV(i)) = 0.0;
+    }
+}
+
+void Mesh::write(void){
+    if(!Utils::writeU("../output/U.txt", U)){
+        std::cout << "Logging: [fem3d] error on writing displacement U" << std::endl;
+        exit(-1);
+    }
+}
+
+// <<< interfaces
+void Mesh::setDBC(void){ // user defined fixed dofs
+    DBCV.resize(0);
+    for(int vI=0; vI<nvet; ++vI){
+        if(abs(vet(vI, 2))<eps){
+            DBCV.conservativeResize(DBCV.size()+3);
+            DBCV(DBCV.size()-3) = vI * 3;
+            DBCV(DBCV.size()-2) = vI * 3 + 1;
+            DBCV(DBCV.size()-1) = vI * 3 + 2;
+        }
+    }
+}
+// >>> interfaces
+
+} // namespace

cpp/src/Mesh.hpp

@@ -0,0 +1,55 @@
+//
+// Mesh.hpp
+//
+// Created by Wei Chen on 9/1/21
+//
+
+#ifndef Mesh_hpp
+#define Mesh_hpp
+
+#include "EigenLibSolver.hpp"
+
+#include <Eigen/Eigen>
+#include <vector>
+#include <set>
+
+namespace fem3d{
+
+class Mesh{
+public: // owned data
+    int nvet, nele, ndof;
+    Eigen::MatrixXd vet; // vertices coordinates
+    Eigen::MatrixXi ele; // vertice index of each tetrahedron
+    Eigen::VectorXd T; // temperatrue of each vertice
+    double E, nu;
+    double alpha, Tref;
+    Eigen::MatrixXd D; // constitutive matrix
+
+public: // owned features
+    double eps = 1e-8;
+    Eigen::VectorXd F; // load of each dof
+    Eigen::VectorXd U; // dofs' displacement to be computed
+    Eigen::VectorXi DBCV; // dofs in DBC
+    std::vector<std::set<int>> vNeighbor; // records all vertices' indices adjacent to each vertice
+    LinSysSolver* linSysSolver;
+
+public: // constructor
+    Mesh(void);
+    ~Mesh(void);
+
+public:
+    void computeFeatures(void);
+    void computeK(void);
+    void computeKe(const Eigen::Matrix<double, 4, 3>& X, double T, 
+                   Eigen::MatrixXd& Ke, Eigen::VectorXd& Fe);
+    void boundaryK(void);
+    void solve(void);
+    void write(void);
+
+public: // interface: set load and DBC
+    void setDBC(void);
+};
+
+} // namespace
+
+#endif // Mesh_hpp

cpp/src/Utils.cpp

@@ -0,0 +1,72 @@
+//
+// Utils.cpp
+//
+// Created by Wei Chen on 9/1/21
+//
+
+#include "Utils.hpp"
+
+namespace fem3d{
+
+bool Utils::readMesh(const std::string& filePath, int& nvet, int& nele,
+              Eigen::MatrixXd& vet, Eigen::MatrixXi& ele){
+    FILE* in = fopen(filePath.c_str(), "r");
+    if(!in){
+        return false;
+    }
+
+    fscanf(in, "%d%d", &nvet, &nele);
+    vet.resize(nvet, 3);
+    ele.resize(nele, 4);
+
+    for(int vI=0; vI<nvet; ++vI){
+        fscanf(in, "%lf%lf%lf", &vet(vI, 0), &vet(vI, 1), &vet(vI, 2));
+    }
+    for(int eleI=0; eleI<nele; ++eleI){
+        fscanf(in, "%d%d%d%d", &ele(eleI, 0), &ele(eleI, 1), &ele(eleI, 2), &ele(eleI, 3));
+    }
+    ele.array() -= 1;
+
+    fclose(in);
+    return true;
+}
+
+bool Utils::readTemp(const std::string& filePath, Eigen::VectorXd& T){
+    FILE* in = fopen(filePath.c_str(), "r");
+    if(!in){
+        return false;
+    }
+
+    for(int vI=0; vI<T.size(); ++vI){
+        fscanf(in, "%lf", &T(vI));
+    }
+
+    fclose(in);
+    return true;
+}
+
+bool Utils::writeU(const std::string& filePath, const Eigen::VectorXd& U){
+    FILE* out = fopen(filePath.c_str(), "w");
+    if(!out){
+        return false;
+    }
+
+    for(int i=0; i<U.size(); ++i){
+        fprintf(out, "%le\n", U(i));
+    }
+
+    fclose(out);
+    return true;
+}
+
+// return sub-matrix(3*3) of H(4*4)
+void Utils::subMatrix(const Eigen::MatrixXd& H, const Eigen::VectorXi& rowInd,
+                      const Eigen::VectorXi& colInd, Eigen::MatrixXd& sub){
+    for(int i=0; i<3; ++i){
+        for(int j=0; j<3; ++j){
+            sub(i, j) = H(rowInd(i), colInd(j));
+        }
+    }
+}
+
+} // namespace

cpp/src/Utils.hpp

@@ -0,0 +1,29 @@
+//
+// Utils.hpp
+//
+// Created by Wei Chen on 9/1/21
+//
+
+#ifndef Utils_hpp
+#define Utils_hpp
+
+#include <Eigen/Eigen>
+#include <iostream>
+#include <fstream>
+
+namespace fem3d{
+
+// a static class implementing some assitant functions
+class Utils{
+public:
+    static bool readMesh(const std::string& filePath, int& nvet, int& nele,
+                         Eigen::MatrixXd& vet, Eigen::MatrixXi& ele);
+    static bool readTemp(const std::string& filePath, Eigen::VectorXd& T);
+    static bool writeU(const std::string& filePath, const Eigen::VectorXd& U);
+    static void subMatrix(const Eigen::MatrixXd& H, const Eigen::VectorXi& rowInd,
+                          const Eigen::VectorXi& colInd, Eigen::MatrixXd& sub);
+};
+
+} // namespace
+
+#endif // Utils_hpp

cpp/src/main.cpp

@@ -0,0 +1,26 @@
+//
+// main.cpp
+//
+// Created by Wei Chen on 8/30/21
+//
+
+#include "Mesh.hpp"
+
+#include <Eigen/Eigen>
+#include <iostream>
+
+int main(){
+    fem3d::Mesh mesh = fem3d::Mesh();
+    mesh.setDBC();
+    std::cout << "Logging: [fem3d] set DBC" << std::endl;
+
+    mesh.computeK();
+    std::cout << "Logging: [fem3d] assembly stiffness matrix and load" << std::endl;
+    mesh.boundaryK();
+
+    mesh.solve();
+    std::cout << "Logging: [fem3d] solve" << std::endl;
+    mesh.write();
+    std::cout << "Logging: [fem3d] write displacement to output/U.txt" << std::endl;
+    return 0;
+}

matlab/fem3d_tet_linear.m

@@ -0,0 +1,152 @@
+%
+% fem3d_tet_linear.m
+% 
+% Created by Wei Chen on 8/4/21
+%
+
+function fem3d_tet_linear()
+nvet = 3168; % nvet: number of vertices, nele: number of elements 
+nele = 16563;
+
+disp('fem simulation start')
+
+eeps = 1e-8;
+% USER-DEFINED MATERIAL PROPERTIES
+E = 73.0;          % Young's modulus of solid material
+nu = 0.17;         % Poisson's ratio 
+alpha = 5.5e-7;
+Tref = 293.15;
+% constitutive matrix
+D = E / (1 + nu) / (1 - 2 * nu) * ...
+[ 1-nu   nu   nu     0          0          0     ;
+    nu 1-nu   nu     0          0          0     ;
+    nu   nu 1-nu     0          0          0     ;
+     0    0    0 (1-2*nu)/2     0          0     ;
+     0    0    0     0      (1-2*nu)/2     0     ;
+     0    0    0     0          0      (1-2*nu)/2];
+
+% read coordinates of vertices, tetrahedron information
+cor = load('input/cor.txt');
+elem = load('input/elem.txt');
+T = load('input/T.txt');
+T = (alpha) * (T - Tref);
+
+% PREPARE FINITE ELEMENT ANALYSIS
+ndof = nvet * 3;
+U = zeros(ndof,1);
+F = zeros(ndof,1);
+
+% USER-DEFINED SUPPORT FIXED DOFs
+fixednids = find(abs(cor(:, 3)) < eeps);
+fixeddofs = [3*fixednids-2; 3*fixednids-1; 3*fixednids];
+freedofs = setdiff(1:ndof,fixeddofs);
+disp('compute DBC')
+
+% assembly the stiffness matrix
+edofMat = kron(elem, 3*ones(1, 3));
+edofMat(:, 1:3:12) = edofMat(:, 1:3:12) - 2;
+edofMat(:, 2:3:12) = edofMat(:, 2:3:12) - 1;
+iK = reshape(kron(edofMat,ones(12,1))',12*12*nele,1);
+jK = reshape(kron(edofMat,ones(1,12))',12*12*nele,1);
+
+sK = zeros(12*12, nele);
+for eleI = 1:nele
+    elenids = [elem(eleI, 1), elem(eleI, 2), elem(eleI, 3), elem(eleI, 4)];
+    eledofs = kron(elenids, 3*ones(1, 3));
+    eledofs(1, 1:3:12) = eledofs(1, 1:3:12) - 2;
+    eledofs(1, 2:3:12) = eledofs(1, 2:3:12) - 1;
+
+    Xe = cor(elenids, :);
+    Te = T(elenids, 1);
+    [KE, Fe] = initKE(D, Xe, Te);
+    sK(:, eleI) = KE(:);
+    F(eledofs, 1) = F(eledofs, 1) + Fe;
+end
+sK = reshape(sK, 12*12*nele, 1);
+
+% FE-ANALYSIS
+K = sparse(iK,jK,sK); K = (K+K')/2;
+% save K1.mat K;
+U(freedofs,:) = K(freedofs,freedofs)\F(freedofs,:);
+disp('solve')
+
+save U.mat U;
+% display result
+cor1 = cor + reshape(U, 3, nvet)';
+figure;
+scatter3(cor1(:, 1), cor1(:, 2), cor1(:, 3));
+
+disp('simulation end')
+
+end
+
+% === GENERATE ELEMENT STIFFNESS MATRIX ===
+% X: 4*3 matrix, which is corridinates of element vertices
+function [KE, Fe] = initKE(D, Xe, Te)
+H = ...
+[ 1.0 Xe(1, 1) Xe(1, 2) Xe(1, 3);
+  1.0 Xe(2, 1) Xe(2, 2) Xe(2, 3);
+  1.0 Xe(3, 1) Xe(3, 2) Xe(3, 3);
+  1.0 Xe(4, 1) Xe(4, 2) Xe(4, 3)];
+V6 = det(H); % 6*V, V is volumn of tetrahedron element
+
+% a1 =  det(H([2, 3, 4], [2, 3, 4]));
+% a2 = -det(H([1, 3, 4], [2, 3, 4]));
+% a3 =  det(H([1, 2, 4], [2, 3, 4]));
+% a4 = -det(H([1, 2, 3], [2, 3, 4]));
+
+b1 = -det(H([2, 3, 4], [1, 3, 4]));
+b2 =  det(H([1, 3, 4], [1, 3, 4]));
+b3 = -det(H([1, 2, 4], [1, 3, 4]));
+b4 =  det(H([1, 2, 3], [1, 3, 4]));
+
+c1 =  det(H([2, 3, 4], [1, 2, 4]));
+c2 = -det(H([1, 3, 4], [1, 2, 4]));
+c3 =  det(H([1, 2, 4], [1, 2, 4]));
+c4 = -det(H([1, 2, 3], [1, 2, 4]));
+
+d1 = -det(H([2, 3, 4], [1, 2, 3]));
+d2 =  det(H([1, 3, 4], [1, 2, 3]));
+d3 = -det(H([1, 2, 4], [1, 2, 3]));
+d4 =  det(H([1, 2, 3], [1, 2, 3]));
+
+B = zeros(6, 12);
+
+B(:, 1:3) = ...
+[ b1  0  0;
+   0 c1  0;
+   0  0 d1;
+  c1 b1  0;
+   0 d1 c1;
+  d1  0 b1];
+
+B(:, 4:6) = ...
+[ b2  0  0;
+   0 c2  0;
+   0  0 d2;
+  c2 b2  0;
+   0 d2 c2;
+  d2  0 b2];
+
+B(:, 7:9) = ...
+[ b3  0  0;
+   0 c3  0;
+   0  0 d3;
+  c3 b3  0;
+   0 d3 c3;
+  d3  0 b3];
+
+B(:, 10:12) = ...
+[ b4  0  0;
+   0 c4  0;
+   0  0 d4;
+  c4 b4  0;
+   0 d4 c4;
+  d4  0 b4];
+
+B = B / V6;
+KE = V6 / 6.0 * (B' * D * B);
+
+strain = mean(Te) * [1; 1; 1; 0; 0; 0];
+Fe = V6 / 6.0 * (B' * D * strain);
+end

python/fem3d_tet_linear.py

@@ -0,0 +1,209 @@
+#
+# fem3d_tet_linear.py
+#
+# Created by Wei Chen on 8/12/21
+#
+
+import numpy as np
+from scipy.sparse import coo_matrix
+from scipy.sparse.linalg import spsolve
+import math
+
+'''
+fem3d_tet_linear: linear fem of 3d tetrahedron model
+@input:
+    nvet: number of vertices
+    nele: number of tetrahedron
+    input/cor.txt: corrdinates of vertices
+    input/elem.txt: indices of tetrahedron elements
+@output:
+    output/u.txt: displacement of vertices
+    output/cor1.txt: result corrdinates of vertices
+'''
+def fem3d_tet_linear():
+    print('fem simulation start')
+
+    eeps = 1e-8
+    # user-defined material properities
+    E = 73.0         # Young's modulus
+    nu = 0.17        # Poisson's ratio
+    alpha = 5.5e-7
+    Tref = 293.15
+    # constitutive matrix
+    D = E / (1. + nu) / (1. - 2. * nu) * np.array(
+        [[1-nu, nu, nu, 0., 0., 0.],
+         [nu, 1-nu, nu, 0., 0., 0.],
+         [nu, nu, 1-nu, 0., 0., 0.],
+         [0., 0., 0., (1.0-2.0*nu)/2.0, 0., 0.],
+         [0., 0., 0., 0., (1.0-2.0*nu)/2.0, 0.],
+         [0., 0., 0., 0., 0., (1.0-2.0*nu)/2.0]], dtype=np.float64)
+    
+    # read coordinates of vertices, tetrahedron information
+    nvet, nele, cor, elem = readTet('input/cube.txt')
+    T = readT(nvet, 'input/cube_T.txt')
+    T = alpha * (T - Tref)
+
+    # prepare finite element analysis
+    ndof = nvet * 3
+    U = np.zeros((ndof), dtype=np.float64)
+    F = np.zeros((ndof), dtype=np.float64)
+
+    # user-defined fixed dofs
+    fixednids = np.asarray(cor[:, 2]==0.0).nonzero()[0]
+    fixeddofs = np.concatenate((fixednids*3, fixednids*3+1, fixednids*3+2), axis=0)
+    freedofs = np.setdiff1d(np.arange(ndof), fixeddofs)
+    print('compute DBC')
+
+    # assembly the stiffness matrix
+    edofMat = np.kron(3*elem, np.ones((1, 3))) + np.kron(np.tile(np.arange(3), 4), np.ones((nele, 1)))
+    edofMat = edofMat.astype(np.int32)
+    iK = np.kron(edofMat, np.ones((12, 1))).flatten()
+    jK = np.kron(edofMat, np.ones((1, 12))).flatten()
+    iK = iK.astype(np.int32)
+    jK = jK.astype(np.int32)
+
+    sK = np.empty((nele, 12*12), dtype=np.float64)
+    for eleI in range(nele):
+        elenids = np.array([elem[eleI, 0], elem[eleI, 1], elem[eleI, 2], elem[eleI, 3]], dtype=np.int32)
+        eledofs = np.kron(3*elenids, np.ones((1, 3))) + np.tile(np.arange(3), 4)
+        eledofs = eledofs.astype(np.int32)
+
+        Xe = cor[elenids, :]
+        Te = T[elenids]
+        KE, Fe = initKE(D, Xe, Te)
+        sK[eleI, :] = KE.flatten('F')
+        F[eledofs] = F[eledofs] + Fe
+    sK = sK.flatten()
+
+    K = coo_matrix((sK, (iK, jK)), shape=(ndof, ndof)).tocsc()
+
+    # FE-analysis
+    print('solve')
+    K = K[freedofs, :][:, freedofs]
+    U[freedofs] = spsolve(K, F[freedofs])
+
+    # write result
+    print('write displacement U.txt and new corrdinates cor1.txt')
+    with open('output/U.txt', 'w') as f:
+        for i in range(ndof):
+            print(U[i], file=f)
+    cor1 = cor + U.reshape(nvet, 3)
+    with open('output/cor1.txt', 'w') as f:
+        for vI in range(nvet):
+            print(cor1[vI, 0], cor1[vI, 1], cor1[vI, 2], file=f)
+
+
+def readTet(filePath):
+    with open(filePath, 'r') as f:
+        line = f.readline()
+        line_split = line.split() 
+        nvet = int(line_split[0])
+        line = f.readline()
+        line_split = line.split()
+        nele = int(line_split[0])
+
+        cor = np.empty((nvet, 3), dtype=np.float64)
+        for vI in range(nvet):
+            line = f.readline()
+            line_split = line.split()
+            cor[vI, 0] = float(line_split[0])
+            cor[vI, 1] = float(line_split[1])
+            cor[vI, 2] = float(line_split[2])
+
+        elem = np.empty((nele, 4), dtype=np.int32)
+        for eleI in range(nele):
+            line = f.readline()
+            line_split = line.split()
+            elem[eleI, 0] = int(line_split[0])
+            elem[eleI, 1] = int(line_split[1])
+            elem[eleI, 2] = int(line_split[2])
+            elem[eleI, 3] = int(line_split[3])
+        elem = elem - 1
+    return nvet, nele, cor, elem
+
+
+def readT(nvet, filePath):
+    T = np.empty((nvet), dtype=np.float64)
+    with open(filePath, 'r') as f:
+        for vI in range(nvet):
+            line = f.readline()
+            line_split = line.split()
+            T[vI] = float(line_split[0])
+    return T
+
+
+# generate element stiffness matrix
+# @input
+# D: constitutive matrix
+# X: 4*3 matrix, which is corridinates of element vertices
+# T: vector of 4 size
+# @output
+# KE: element stiffness matrix
+def initKE(D, X, T):
+    H = np.array([[1., X[0, 0], X[0, 1], X[0, 2]],
+                  [1., X[1, 0], X[1, 1], X[1, 2]],
+                  [1., X[2, 0], X[2, 1], X[2, 2]],
+                  [1., X[3, 0], X[3, 1], X[3, 2]]])
+    V6 = abs(np.linalg.det(H))
+
+    a1 =  np.linalg.det(H[[1, 2, 3], :][:, [1, 2, 3]])    
+    a2 = -np.linalg.det(H[[0, 2, 3], :][:, [1, 2, 3]])    
+    a3 =  np.linalg.det(H[[0, 1, 3], :][:, [1, 2, 3]])    
+    a4 = -np.linalg.det(H[[0, 1, 2], :][:, [1, 2, 3]])    
+
+    b1 = -np.linalg.det(H[[1, 2, 3], :][:, [0, 2, 3]])    
+    b2 =  np.linalg.det(H[[0, 2, 3], :][:, [0, 2, 3]])    
+    b3 = -np.linalg.det(H[[0, 1, 3], :][:, [0, 2, 3]])    
+    b4 =  np.linalg.det(H[[0, 1, 2], :][:, [0, 2, 3]])    
+
+    c1 =  np.linalg.det(H[[1, 2, 3], :][:, [0, 1, 3]])    
+    c2 = -np.linalg.det(H[[0, 2, 3], :][:, [0, 1, 3]])    
+    c3 =  np.linalg.det(H[[0, 1, 3], :][:, [0, 1, 3]])    
+    c4 = -np.linalg.det(H[[0, 1, 2], :][:, [0, 1, 3]])    
+
+    d1 = -np.linalg.det(H[[1, 2, 3], :][:, [0, 1, 2]])    
+    d2 =  np.linalg.det(H[[0, 2, 3], :][:, [0, 1, 2]])    
+    d3 = -np.linalg.det(H[[0, 1, 3], :][:, [0, 1, 2]])    
+    d4 =  np.linalg.det(H[[0, 1, 2], :][:, [0, 1, 2]])    
+
+    B = np.empty((6, 12), dtype=np.float64)
+
+    B[:, 0:3] = np.array([[b1, 0., 0.],
+                          [0., c1, 0.],
+                          [0., 0., d1],
+                          [c1, b1, 0.],
+                          [0., d1, c1],
+                          [d1, 0., b1]])
+
+    B[:, 3:6] = np.array([[b2, 0., 0.],
+                          [0., c2, 0.],
+                          [0., 0., d2],
+                          [c2, b2, 0.],
+                          [0., d2, c2],
+                          [d2, 0., b2]])
+
+    B[:, 6:9] = np.array([[b3, 0., 0.],
+                          [0., c3, 0.],
+                          [0., 0., d3],
+                          [c3, b3, 0.],
+                          [0., d3, c3],
+                          [d3, 0., b3]])
+
+    B[:, 9:12] = np.array([[b4, 0., 0.],
+                           [0., c4, 0.],
+                           [0., 0., d4],
+                           [c4, b4, 0.],
+                           [0., d4, c4],
+                           [d4, 0., b4]])
+    
+    B = B / V6
+    KE = V6 / 6.0 * (B.T @ D @ B)
+    
+    strain = np.mean(T) * np.array([1., 1., 1., 0., 0., 0.])
+    Fe = V6 / 6.0 * (B.T @ D @ strain)
+
+    return KE, Fe
+
+
+if __name__ == '__main__':
+    fem3d_tet_linear()