<repo_name>zhengzhirun/Finite-Element-Methods<file_sep>/Possion方程--C++(CVT自适应)/Possion.h
/*************************************************************************
	> File Name: Possion.h
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年08月19日 星期日 22时04分12秒
 ************************************************************************/

#ifndef _POSSION_H
#define _POSSION_H

#include "Matrix.h"

double Possion(const std::string, const std::string, std::vector<double>& ,int);    // 求解Possion方程
void nodesf2dat(const std::string, std::vector<Point>&, std::vector<int>&);    // 读取nodes.dat文档
void trigsf2dat(const std::string, std::vector<std::vector<int>>&);     // 读取trigs.dat文档
void writeTriangle(std::vector<Point>&, std::vector<std::vector<int>>&, std::vector<int>&, 
        const std::string);

#endif
<file_sep>/Possion方程--C++(CVT自适应)/Matrix.h
// 设置类的保护符
#ifndef MATRIX_H
#define MATRIX_H

#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <fstream>
#include <Eigen/Dense>  // 线性代数库
#include <Eigen/Sparse> // 稀疏存储矩阵库,以及稀疏存储的线性方程组求解库
#include <Eigen/IterativeLinearSolvers> // 迭代法求解线性方程组的求解库
#include <unsupported/Eigen/IterativeSolvers>

class Point
{
	friend std::ostream& operator<< (std::ostream&, const Point&);
	friend std::istream& operator>> (std::istream&, Point&);
	friend Point operator+ (const Point&, const Point&);
	friend Point operator- (const Point&, const Point&);
	friend Point midpoint(const Point&, const Point&);
	friend double distance(const Point&, const Point&);
    // 求三角形的质心,可以求加权质心
	friend Point barycenter(const std::vector<Point>&, const double *);
    friend double AREA(const Point&, const Point&, const Point&);

	private:
		double x[2];   
	
	public:
		Point();  
		Point(const double*);  
 		Point(const Point&);
 		Point(const double&, const double&);  

 		Point& operator= (const Point&);
 		double& operator[] (int);
 		const double& operator[] (int) const;
 		double length() const;
 		Point& operator+= (const Point&);		
 		Point& operator-= (const Point&);
 		Point& operator*= (const double&);
 		Point& operator/= (const double&);
};

class TmpEle
{
	friend std::ostream& operator<< (std::ostream&, const TmpEle&);

	private:
		std::vector<Point> Pnt;
		std::vector<Point> GaussPnt;
		std::vector<double> GaussWeight;

	public:
		// 构造函数
		TmpEle();
		TmpEle(int);
		void buildTE(int);

		Point getPnt(int);
		double getArea();
		std::vector<Point> getGaussPnt();
		std::vector<double> getGaussWeight();

		// 标准单元到普通单元
		Point Local_to_Global(const Point&, const std::vector<Point> &) const;
		// 普通单元到标准单元
		Point Global_to_Local(const Point&, const std::vector<Point> &) const;
		// 标准单元到普通单元
		std::vector<Point> Local_to_Global(const std::vector<Point> &,const std::vector<Point> & ) const;
		// 普通单元到标准单元
		std::vector<Point> Global_to_Local(const std::vector<Point>&, const std::vector<Point>&) const;
		// 从标准单元到普通单元的雅克比矩阵
		double Local_to_Global_jacobian(const Point&, const std::vector<Point>&) const;
		// 从普通单元到标准单元的雅克比矩阵
		double Global_to_Local_jacobian(const Point&, const std::vector<Point>&) const;
		// 从标准单元到普通单元的雅克比矩阵
		std::vector<double> Local_to_Global_jacobian(const std::vector<Point>&, const std::vector<Point>&) const;
		// 从普通单元到标准单元的雅克比矩阵
		std::vector<double> Global_to_Local_jacobian(const std::vector<Point>&, const std::vector<Point>&) const; 
};

class Mesh
{
	friend std::ostream& operator<< (std::ostream&, const Mesh&);

	private:
		std::vector<Point> Pnt;   // 三角单元的点
		std::vector<int> BndPnt;  // 三角单元的边界点
		std::vector<std::vector<int>> Ele;  // 三角单元的边

	public:
		Mesh() = default;  // 默认构造函数

		void readData(const std::string&);
		int getEleVtx(int , int);
		std::vector<int> getEleVtx(int);
		Point getPnt(int);
		std::vector<Point> getPnt(std::vector<int>&);
		std::vector<int> getBndPnt();
		int n_element(); // 得到三角单元的数量
		int n_point();  // 得到点的个数
		int n_boundaryPoint(); // 得到边界点的个数
        void neighbor(std::vector<std::vector<int>>&);  // 单元邻居索引
};

class PDE
{	
	public:
		// 偏微分方程的边界条件
		double u_boundary(const Point&);

		// 偏微分方程的右端项
		double f(const Point&);

		// 偏微分方程解析解
		double u_exact(const Point&);
};


class Matrix
{
	public:
		// 数据结构
		PDE pde;	
		Mesh mesh;     
		TmpEle tmpEle;  
		std::vector<double> u_h;

		Matrix() = default;
		Matrix(const std::string&, int);  
		// gauss-seidel迭代
		void GaussSeidel(const std::vector<std::vector<double>>&,std::vector<double>&, 
                   const std::vector<double>&);
        // 利用Eigen库求解稀疏代数方程组
        void SolveLinearASparse(const Eigen::SparseMatrix<double>&,std::vector<double>&, 
                                const std::vector<double>&, int);
		// 基函数的构造
		void basisValue(const Point&, const std::vector<Point>&, std::vector<double>&);
		void basisValue(const std::vector<Point>&, const std::vector<Point>&, std::vector<std::vector<double>>&);
		// 梯度算子的计算
		void basisGrad(const Point&, const std::vector<Point>&, std::vector<std::vector<double>>&);
		void basisGrad(const std::vector<Point>&, const std::vector<Point>&, 
                 std::vector<std::vector<std::vector<double>>>&);
		// 单元刚度矩阵
		void buildEleMatrix_RHS(int, std::vector<std::vector<double>>&, std::vector<double>&);	    
        // 总刚度矩阵(稀疏存储)
        void buildMatrix_RHS_Sparse(Eigen::SparseMatrix<double>&, std::vector<double>&);        
        // 稀疏矩阵的边界条件的处理
        void DealBoundary_Sparse(Eigen::SparseMatrix<double>&, const std::vector<double>&, std::vector<double>&);
		// 计算L2误差
		double ComputerL2Error(const std::vector<double> &);        
        // 计算剖分上每一个单元上的梯度
        void GradBasis(std::vector<std::vector<std::vector<double>>>&);
        // 计算u_h的梯度
        void Gradu(const std::vector<std::vector<std::vector<double>>>&, 
                   const std::vector<double>&, std::vector<std::vector<double>>&);
        // 三角单元上每一条边的向量
        void eleVector(std::vector<std::vector<std::vector<double>>>&);
        // 三角单元上每一条边的法向量
        void normalVector(const std::vector<std::vector<std::vector<double>>>&,
                      std::vector<std::vector<std::vector<double>>>&);
        // 三角单元的直径
        double diameter(const std::vector<Point>&);
        // residual type(误差指示子)
        void Estimateresidual(std::vector<double>&,std::vector<double>&); 
};


#endif
<file_sep>/Possion方程--C++(CVT自适应)/density.h
/*************************************************************************
	> File Name: density.h
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年08月20日 星期一 22时37分18秒
 ************************************************************************/

#ifndef _DENSITY_H
#define _DENSITY_H

#include<fstream>
#include<vector>
#include<string>
#include "Matrix.h"

void density(const std::vector<Point>&, const std::vector<std::vector<int>>&, const std::vector<double>&);


#endif
<file_sep>/CVT网格控制-绘制误差/CVT_Error.cpp
/*************************************************************************
	> File Name: CVT_Error.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年08月23日 星期四 17时13分50秒
 ************************************************************************/
#include "Possion.h"
int main()
{
    // 网格初始化
    system("./mesh_init -f square.poly -q 1 -a 0.01");
    // 把初始网格信息写入文档中
    std::vector<Point> node;
    std::vector<std::vector<int>> ele;
    std::vector<int> boundary_node;
    writeTriangle(node,ele,boundary_node,"Trimesh_0.txt");  // 初始网格信息
    std::vector<double> Error;  // Error[0]为L2范误差,Error[1]为H1范误差
    std::vector<std::vector<double>> Errors; // 存储每一次的误差
    std::vector<double> eta;
    Error = Possion("Trimesh_0.txt","Results.txt",eta,7);
    // 记录点的个数
    std::vector<int> nodes; 
    system("./dens_info");
    Errors.push_back(Error);
    nodes.push_back(node.size());
    system("./mesh_opti -b 1");

    // 文件处理
    for (int i = 0; i != 200; ++i){
        std::string filename = "Trimesh_";
        filename = filename + std::to_string(i+1) + ".txt";
        writeTriangle(node,ele,boundary_node,filename,i);
        nodes.push_back(node.size());
        std::cout << "i = " << i << std::endl;
        Error = Possion(filename,"Results.txt",eta,7);
        Errors.push_back(Error);
    }
    // 输出文档
    std::ofstream out("Error.txt",std::ofstream::out);    
    for (decltype(Errors.size()) i = 0; i != Errors.size(); ++i)
        out << i+1 << "\t" << Errors[i][0] << "\t\t" << Errors[i][1] << "\t\t" << nodes[i] << "\n";
    out.close();
    
    return 0;
}
<file_sep>/CVT网格控制-绘制误差/Makefile
all: CVT_Error

CVT_Error: CVT_Error.o Matrix.o Possion.o
	g++ -std=c++0x -o CVT_Error CVT_Error.o Possion.o Matrix.o 

CVT_Error.o: CVT_Error.cpp Possion.h 
	g++ -std=c++0x -c CVT_Error.cpp

Possion.o: Possion.cpp Possion.h
	g++ -std=c++0x -c Possion.cpp

Matrix.o: Matrix.cpp Matrix.h
	g++ -std=c++0x -c Matrix.cpp
<file_sep>/Heat方程-C++隐格式-非零边界/Mesh.h
/*************************************************************************
	> File Name: Mesh.h
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时35分13秒
 ************************************************************************/

#ifndef MESH_H
#define MESH_H

#include "Point.h"
#include <fstream>
class Mesh
{
friend std::ostream& operator<< (std::ostream&, const Mesh&);

private:
	std::vector<Point> Pnt;   // 三角单元的点
	std::vector<int> BndPnt;  // 三角单元的边界点
	std::vector<std::vector<int>> Ele;  // 三角单元的边

public:
	Mesh() = default;  // 默认构造函数
	void readData(const std::string&);
	int getEleVtx(int , int);
	std::vector<int> getEleVtx(int);
	Point getPnt(int);
	std::vector<Point> getPnt(std::vector<int>&);
	std::vector<int> getBndPnt();
	int n_element(); // 得到三角单元的数量
	int n_point();  // 得到点的个数
	int n_boundaryPoint(); // 得到边界点的个数
};

#endif
<file_sep>/Heat方程-C++隐格式/TmpEle.cpp
/*************************************************************************
	> File Name: TmpEle.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时10分27秒
 ************************************************************************/

#include "TmpEle.h"
// 非成员函数
std::ostream& operator<< (std::ostream& os, const TmpEle& te)
{	
	for (auto i : te.Pnt)
		os << i << '\n';
	for (auto i : te.GaussPnt)
		os << i << '\n';
	for (auto i : te.GaussWeight)
		os << i << '\n';
	return os;
}

// 构造函数
TmpEle::TmpEle()
{
	Pnt.resize(3);
	Pnt[0][0]=0.0;
	Pnt[0][1]=0.0; 
	Pnt[1][0]=1.0;
	Pnt[1][1]=0.0; 
	Pnt[2][0]=0.0;
	Pnt[2][1]=1.0; 
	
	GaussPnt.resize(7);
	GaussWeight.resize(7);
	GaussPnt[0][0]=0.333333333333;
	GaussPnt[1][0]=0.470142064105;
	GaussPnt[2][0]=0.470142064105;
	GaussPnt[3][0]=0.059715871790;
	GaussPnt[4][0]=0.101286507323;
	GaussPnt[5][0]=0.101286507323;
	GaussPnt[6][0]=0.797426985353;
	GaussPnt[0][1]=0.333333333333;
	GaussPnt[1][1]=0.470142064105;
	GaussPnt[2][1]=0.059715871790;
	GaussPnt[3][1]=0.470142064105;
	GaussPnt[4][1]=0.101286507323;
	GaussPnt[5][1]=0.797426985353;
	GaussPnt[6][1]=0.101286507323;
	GaussWeight[0]=0.225000000000;
	GaussWeight[1]=0.132394152788;
	GaussWeight[2]=0.132394152789;
	GaussWeight[3]=0.132394152788;
	GaussWeight[4]=0.125939180545;
	GaussWeight[5]=0.125939180545; 
	GaussWeight[6]=0.125939180545; 
}

TmpEle::TmpEle(int i)
{
  Pnt.resize(3);
  Pnt[0][0]=0.0;
  Pnt[0][1]=0.0; 
  Pnt[1][0]=1.0;
  Pnt[1][1]=0.0; 
  Pnt[2][0]=0.0;
  Pnt[2][1]=1.0; 
  switch(i){
     case 1:
        GaussPnt.resize(1);
		GaussWeight.resize(1);
		GaussPnt[0][0]=0.3333333333333333;
		GaussPnt[0][1]=0.3333333333333333;
		GaussWeight[0]=1.0;
        break;
     case 2:
        GaussPnt.resize(3);
		GaussWeight.resize(3);
		GaussPnt[0][0]=0.166666666667;
		GaussPnt[1][0]=0.166666666667;
		GaussPnt[2][0]=0.666666666667;
		GaussPnt[0][1]=0.166666666667;
		GaussPnt[1][1]=0.666666666667;
		GaussPnt[2][1]=0.166666666667;
		GaussWeight[0]=0.333333333333;
		GaussWeight[1]=0.333333333334;
		GaussWeight[2]=0.333333333333;
        break;
     case 3:
        GaussPnt.resize(4);
		GaussWeight.resize(4);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.200000000000;
		GaussPnt[2][0]=0.200000000000;
		GaussPnt[3][0]=0.600000000000;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.200000000000;
		GaussPnt[2][1]=0.600000000000;
		GaussPnt[3][1]=0.200000000000;
		GaussWeight[0]=-0.562500000000;
		GaussWeight[1]=0.520833333333;
		GaussWeight[2]=0.520833333334;
		GaussWeight[3]=0.520833333333;
        break;
     case 4:
        GaussPnt.resize(6);
		GaussWeight.resize(6);
		GaussPnt[0][0]=0.445948490916;
		GaussPnt[1][0]=0.445948490916;
		GaussPnt[2][0]=0.108103018168;
		GaussPnt[3][0]=0.091576213510;
		GaussPnt[4][0]=0.091576213510;
		GaussPnt[5][0]=0.816847572980;
		GaussPnt[0][1]=0.445948490916;
		GaussPnt[1][1]=0.108103018168;
		GaussPnt[2][1]=0.445948490916;
		GaussPnt[3][1]=0.091576213510;
		GaussPnt[4][1]=0.816847572980;
		GaussPnt[5][1]=0.091576213510;
		GaussWeight[0]=0.223381589678;
		GaussWeight[1]=0.223381589678;
		GaussWeight[2]=0.223381589678;
		GaussWeight[3]=0.109951743655;
		GaussWeight[4]=0.109951743656;
		GaussWeight[5]=0.109951743655;
        break;
     case 5:
        GaussPnt.resize(7);
		GaussWeight.resize(7);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.470142064105;
		GaussPnt[2][0]=0.470142064105;
		GaussPnt[3][0]=0.059715871790;
		GaussPnt[4][0]=0.101286507323;
		GaussPnt[5][0]=0.101286507323;
		GaussPnt[6][0]=0.797426985353;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.470142064105;
		GaussPnt[2][1]=0.059715871790;
		GaussPnt[3][1]=0.470142064105;
		GaussPnt[4][1]=0.101286507323;
		GaussPnt[5][1]=0.797426985353;
		GaussPnt[6][1]=0.101286507323;
		GaussWeight[0]=0.225000000000;
		GaussWeight[1]=0.132394152788;
		GaussWeight[2]=0.132394152789;
		GaussWeight[3]=0.132394152788;
		GaussWeight[4]=0.125939180545;
		GaussWeight[5]=0.125939180545; 
		GaussWeight[6]=0.125939180545; 
        break;
     case 6:
        GaussPnt.resize(12);
		GaussWeight.resize(12);
		GaussPnt[0][0]=0.249286745171;
		GaussPnt[1][0]=0.249286745171;
		GaussPnt[2][0]=0.501426509658;
		GaussPnt[3][0]=0.063089014492;
		GaussPnt[4][0]=0.063089014492;
		GaussPnt[5][0]=0.873821971017;
		GaussPnt[6][0]=0.310352451034;
		GaussPnt[7][0]=0.636502499121;
		GaussPnt[8][0]=0.053145049845;
		GaussPnt[9][0]=0.636502499121;
		GaussPnt[10][0]=0.310352451034;
		GaussPnt[11][0]=0.053145049845;
		GaussPnt[0][1]=0.249286745171;
		GaussPnt[1][1]=0.501426509658;
		GaussPnt[2][1]=0.249286745171;
		GaussPnt[3][1]=0.063089014492;
		GaussPnt[4][1]=0.873821971017;
		GaussPnt[5][1]=0.063089014492;
		GaussPnt[6][1]=0.636502499121;
		GaussPnt[7][1]=0.053145049845;
		GaussPnt[8][1]=0.310352451034;
		GaussPnt[9][1]=0.310352451034;
		GaussPnt[10][1]=0.053145049845;
		GaussPnt[11][1]=0.636502499121;
		GaussWeight[0]=0.116786275726;
		GaussWeight[1]=0.116786275726;
		GaussWeight[2]=0.116786275726;
		GaussWeight[3]=0.050844906370;
		GaussWeight[4]=0.050844906370;
		GaussWeight[5]=0.050844906370; 
		GaussWeight[6]=0.082851075618;
		GaussWeight[7]=0.082851075618; 
		GaussWeight[8]=0.082851075618; 
		GaussWeight[9]=0.082851075618; 
		GaussWeight[10]=0.082851075618; 
		GaussWeight[11]=0.082851075618; 
        break;
     case 7:
        GaussPnt.resize(13);
		GaussWeight.resize(13);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.260345966079;
		GaussPnt[2][0]=0.260345966079;
		GaussPnt[3][0]=0.479308067842;
		GaussPnt[4][0]=0.065130102902;
		GaussPnt[5][0]=0.065130102902;
		GaussPnt[6][0]=0.869739794196;
		GaussPnt[7][0]=0.312865496005;
		GaussPnt[8][0]=0.638444188570;
		GaussPnt[9][0]=0.048690315425;
		GaussPnt[10][0]=0.638444188570;
		GaussPnt[11][0]=0.312865496005;
		GaussPnt[12][0]=0.048690315425;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.260345966079;
		GaussPnt[2][1]=0.479308067842;
		GaussPnt[3][1]=0.260345966079;
		GaussPnt[4][1]=0.065130102902;
		GaussPnt[5][1]=0.869739794196;
		GaussPnt[6][1]=0.065130102902;
		GaussPnt[7][1]=0.638444188570;
		GaussPnt[8][1]=0.048690315425;
		GaussPnt[9][1]=0.312865496005;
		GaussPnt[10][1]=0.312865496005;
		GaussPnt[11][1]=0.048690315425;
		GaussPnt[12][1]=0.638444188570;
		GaussWeight[0]=-0.149570044468;
		GaussWeight[1]=0.175615257433;
		GaussWeight[2]=0.175615257433;
		GaussWeight[3]=0.175615257433;
		GaussWeight[4]=0.053347235609;
		GaussWeight[5]=0.053347235609; 
		GaussWeight[6]=0.053347235609;
		GaussWeight[7]=0.077113760890; 
		GaussWeight[8]=0.077113760890; 
		GaussWeight[9]=0.077113760890; 
		GaussWeight[10]=0.077113760890; 
		GaussWeight[11]=0.077113760890; 
		GaussWeight[12]=0.077113760890;  
        break;
     case 10:
        GaussPnt.resize(25);
		GaussWeight.resize(25);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.485577633384;
		GaussPnt[2][0]=0.485577633384;
		GaussPnt[3][0]=0.028844733233;
		GaussPnt[4][0]=0.109481575485;
		GaussPnt[5][0]=0.109481575485;
		GaussPnt[6][0]=0.781036849030;
		GaussPnt[7][0]=0.307939838764;
		GaussPnt[8][0]=0.550352941821;
		GaussPnt[9][0]=0.141707219415;
		GaussPnt[10][0]=0.550352941821;
		GaussPnt[11][0]=0.307939838764;
		GaussPnt[12][0]=0.141707219415;
		GaussPnt[13][0]=0.246672560640;
		GaussPnt[14][0]=0.728323904597;
		GaussPnt[15][0]=0.025003534763;
		GaussPnt[16][0]=0.728323904597;
		GaussPnt[17][0]=0.246672560640;
		GaussPnt[18][0]=0.025003534763;
		GaussPnt[19][0]=0.066803251012;
		GaussPnt[20][0]=0.923655933587;
		GaussPnt[21][0]=0.009540815400;
		GaussPnt[22][0]=0.923655933587;
		GaussPnt[23][0]=0.066803251012;
		GaussPnt[24][0]=0.009540815400;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.485577633384;
		GaussPnt[2][1]=0.028844733233;
		GaussPnt[3][1]=0.485577633384;
		GaussPnt[4][1]=0.109481575485;
		GaussPnt[5][1]=0.781036849030;
		GaussPnt[6][1]=0.109481575485;
		GaussPnt[7][1]=0.550352941821;
		GaussPnt[8][1]=0.141707219415;
		GaussPnt[9][1]=0.307939838764;
		GaussPnt[10][1]=0.307939838764;
		GaussPnt[11][1]=0.141707219415;
		GaussPnt[12][1]=0.550352941821;
		GaussPnt[13][1]=0.728323904597;
		GaussPnt[14][1]=0.025003534763;
		GaussPnt[15][1]=0.246672560640;
		GaussPnt[16][1]=0.246672560640;
		GaussPnt[17][1]=0.025003534763;
		GaussPnt[18][1]=0.728323904597;
		GaussPnt[19][1]=0.923655933587;
		GaussPnt[20][1]=0.009540815400;
		GaussPnt[21][1]=0.066803251012;
		GaussPnt[22][1]=0.066803251012;
		GaussPnt[23][1]=0.009540815400;
		GaussPnt[24][1]=0.923655933587;
		GaussWeight[0]=0.090817990383;
		GaussWeight[1]=0.036725957756;
		GaussWeight[2]=0.036725957756;
		GaussWeight[3]=0.036725957756;
		GaussWeight[4]=0.045321059436;
		GaussWeight[5]=0.045321059436; 
		GaussWeight[6]=0.045321059436;
		GaussWeight[7]=0.072757916845; 
		GaussWeight[8]=0.072757916845; 
		GaussWeight[9]=0.072757916845; 
		GaussWeight[10]=0.072757916845; 
		GaussWeight[11]=0.072757916845; 
		GaussWeight[12]=0.072757916845; 
		GaussWeight[13]=0.028327242531; 
		GaussWeight[14]=0.028327242531; 
		GaussWeight[15]=0.028327242531; 
		GaussWeight[16]=0.028327242531; 
		GaussWeight[17]=0.028327242531; 
		GaussWeight[18]=0.028327242531; 
		GaussWeight[19]=0.009421666964; 
		GaussWeight[20]=0.009421666964; 
		GaussWeight[21]=0.009421666964; 
		GaussWeight[22]=0.009421666964; 
		GaussWeight[23]=0.009421666964; 
		GaussWeight[24]=0.009421666964; 
        break;
     case 12:
        GaussPnt.resize(33);
		GaussWeight.resize(33);
		GaussPnt[0][0]=0.488217389774;
		GaussPnt[1][0]=0.488217389774;
		GaussPnt[2][0]=0.023565220452;
		GaussPnt[3][0]=0.439724392294;
		GaussPnt[4][0]=0.439724392294;
		GaussPnt[5][0]=0.120551215411;
		GaussPnt[6][0]=0.271210385012;
		GaussPnt[7][0]=0.271210385012;
		GaussPnt[8][0]=0.457579229976;
		GaussPnt[9][0]=0.127576145542;
		GaussPnt[10][0]=0.127576145542;
		GaussPnt[11][0]=0.744847708917;
		GaussPnt[12][0]=0.021317350453;
		GaussPnt[13][0]=0.021317350453;
		GaussPnt[14][0]=0.957365299094;
		GaussPnt[15][0]=0.275713269686;
		GaussPnt[16][0]=0.608943235780;
		GaussPnt[17][0]=0.115343494535;
		GaussPnt[18][0]=0.608943235780;
		GaussPnt[19][0]=0.275713269686;
		GaussPnt[20][0]=0.115343494535;
		GaussPnt[21][0]=0.281325580990;
		GaussPnt[22][0]=0.695836086788;
		GaussPnt[23][0]=0.022838332222;
		GaussPnt[24][0]=0.695836086788;
		GaussPnt[25][0]=0.281325580990;
		GaussPnt[26][0]=0.022838332222;
		GaussPnt[27][0]=0.116251915908;
		GaussPnt[28][0]=0.858014033544;
		GaussPnt[29][0]=0.025734050548;
		GaussPnt[30][0]=0.858014033544;
		GaussPnt[31][0]=0.116251915908;
		GaussPnt[32][0]=0.025734050548;
		GaussPnt[0][1]=0.488217389774;
		GaussPnt[1][1]=0.023565220452;
		GaussPnt[2][1]=0.488217389774;
		GaussPnt[3][1]=0.439724392294;
		GaussPnt[4][1]=0.120551215411;
		GaussPnt[5][1]=0.439724392294;
		GaussPnt[6][1]=0.271210385012;
		GaussPnt[7][1]=0.457579229976;
		GaussPnt[8][1]=0.271210385012;
		GaussPnt[9][1]=0.127576145542;
		GaussPnt[10][1]=0.744847708917;
		GaussPnt[11][1]=0.127576145542;
		GaussPnt[12][1]=0.021317350453;
		GaussPnt[13][1]=0.957365299094;
		GaussPnt[14][1]=0.021317350453;
		GaussPnt[15][1]=0.608943235780;
		GaussPnt[16][1]=0.115343494535;
		GaussPnt[17][1]=0.275713269686;
		GaussPnt[18][1]=0.275713269686;
		GaussPnt[19][1]=0.115343494535;
		GaussPnt[20][1]=0.608943235780;
		GaussPnt[21][1]=0.695836086788;
		GaussPnt[22][1]=0.022838332222;
		GaussPnt[23][1]=0.281325580990;
		GaussPnt[24][1]=0.281325580990;
		GaussPnt[25][1]=0.022838332222;
		GaussPnt[26][1]=0.695836086788;
		GaussPnt[27][1]=0.858014033544;
		GaussPnt[28][1]=0.025734050548;
		GaussPnt[29][1]=0.116251915908;
		GaussPnt[30][1]=0.116251915908;
		GaussPnt[31][1]=0.025734050548;
		GaussPnt[32][1]=0.858014033544;
		GaussWeight[0]=0.025731066440;
		GaussWeight[1]=0.025731066440;
		GaussWeight[2]=0.025731066440;
		GaussWeight[3]=0.043692544538;
		GaussWeight[4]=0.043692544538;
		GaussWeight[5]=0.043692544538; 
		GaussWeight[6]=0.062858224218;
		GaussWeight[7]=0.062858224218; 
		GaussWeight[8]=0.062858224218; 
		GaussWeight[9]=0.034796112931; 
		GaussWeight[10]=0.034796112931; 
		GaussWeight[11]=0.034796112931; 
		GaussWeight[12]=0.006166261052; 
		GaussWeight[13]=0.006166261052; 
		GaussWeight[14]=0.006166261052; 
		GaussWeight[15]=0.040371557766; 
		GaussWeight[16]=0.040371557766; 
		GaussWeight[17]=0.040371557766; 
		GaussWeight[18]=0.040371557766; 
		GaussWeight[19]=0.040371557766; 
		GaussWeight[20]=0.040371557766; 
		GaussWeight[21]=0.022356773202; 
		GaussWeight[22]=0.022356773202; 
		GaussWeight[23]=0.022356773202; 
		GaussWeight[24]=0.022356773202; 
		GaussWeight[25]=0.022356773202; 
		GaussWeight[26]=0.022356773202; 
		GaussWeight[27]=0.017316231109; 
		GaussWeight[28]=0.017316231109; 
		GaussWeight[29]=0.017316231109; 
		GaussWeight[30]=0.017316231109; 
		GaussWeight[31]=0.017316231109; 
		GaussWeight[32]=0.017316231109; 
        break;
     default:
        std::cout<<"There are no Gauss quadrature points with accuracy of "<<i<<"-th order"<<std::endl;
  }
}

// 成员函数

void TmpEle::buildTE(int i)
{
	TmpEle tmpEle(i);
	for (auto x : tmpEle.Pnt)
		Pnt.push_back(x);
	for (auto x : tmpEle.GaussPnt)
		GaussPnt.push_back(x);
	for (auto x : tmpEle.GaussWeight)
		GaussWeight.push_back(x);
}

Point TmpEle::getPnt(int i)
{
	return Pnt[i];
}

std::vector<Point> TmpEle::getGaussPnt()
{
	return GaussPnt;
}

std::vector<double> TmpEle::getGaussWeight()
{
	return GaussWeight;
}

double TmpEle::getArea()
{
	return 0.5*((Pnt[1][0] - Pnt[0][0]) * (Pnt[2][1] - Pnt[0][1])
				- (Pnt[1][1] - Pnt[0][1]) * (Pnt[2][0] - Pnt[0][0]));
}

// 标准单元到普通单元
Point TmpEle::Local_to_Global(const Point & lp, const std::vector<Point> & gv) const
{
	Point gp;
	double lambda[3];
	double area = AREA(Pnt[0], Pnt[1], Pnt[2]);
	lambda[0] = AREA(lp, Pnt[1], Pnt[2]) / area;
	lambda[1] = AREA(lp, Pnt[2], Pnt[0]) / area;
	lambda[2] = AREA(lp, Pnt[0], Pnt[1]) / area;
	gp[0] = lambda[0] * gv[0][0] + lambda[1] * gv[1][0] + lambda[2] * gv[2][0];
	gp[1] = lambda[0] * gv[0][1] + lambda[1] * gv[1][1] + lambda[2] * gv[2][1];
	return gp;
}

// 普通单元到标准单元
Point TmpEle::Global_to_Local(const Point& gp, const std::vector<Point>& gv) const
{
	Point lp;
	double area = AREA(gv[0],gv[1],gv[2]);

	lp[0] = ((gv[2][1] - gv[0][1]) * (gp[0] - gv[0][0]) - 
		(gv[2][0] - gv[0][0]) * (gp[1] - gv[0][1])) / (2 * area);
	lp[1] = (-(gv[1][1] - gv[0][1]) * (gp[0] - gv[0][0]) +
		(gv[1][0] - gv[0][0]) * (gp[1] - gv[0][1])) / (2 * area);

	return lp;
}

// 标准单元到普通单元
std::vector<Point> TmpEle::Local_to_Global(const std::vector<Point>& lp, const std::vector<Point> & gv) const
{
	std::vector<Point> gp(lp.size());
	double area = AREA(Pnt[0], Pnt[1], Pnt[2]);
	for(int i=0; i<gp.size(); i++){
		double lambda[3];
		lambda[0] = AREA(lp[i], Pnt[1], Pnt[2]) / area;
		lambda[1] = AREA(lp[i], Pnt[2], Pnt[0]) / area;
		lambda[2] = AREA(lp[i], Pnt[0], Pnt[1]) / area;
		gp[i][0] = lambda[0] * gv[0][0] + lambda[1] * gv[1][0] + lambda[2] * gv[2][0];
		gp[i][1] = lambda[0] * gv[0][1] + lambda[1] * gv[1][1] + lambda[2] * gv[2][1];
	}
	return gp;
}

// 普通单元到标准单元
std::vector<Point> TmpEle::Global_to_Local(const std::vector<Point>& gp, 
	const std::vector<Point>& gv) const
{
	std::vector<Point> lp;
	double area = AREA(gv[0],gv[1],gv[2]);

	for (auto point : gp){
		Point p;
		p[0] = ((gv[2][1] - gv[0][1]) * (point[0] - gv[0][0]) - 
				(gv[2][0] - gv[0][0]) * (point[1] - gv[0][1])) / (2 * area);
		p[1] = (-(gv[1][1] - gv[0][1]) * (point[0] - gv[0][0]) + 
				(gv[1][0] - gv[0][0]) * (point[1] - gv[0][1])) / (2 * area);
		lp.push_back(p);
	}
	return lp;
}

// 从标准单元到普通单元的雅克比矩阵
double TmpEle::Local_to_Global_jacobian(const Point& lp, const std::vector<Point>& gv) const
{
	// lp只是提供维度
	double jacobian =  AREA(gv[0],gv[1],gv[2]) / AREA(Pnt[0], Pnt[1], Pnt[2]);
	return jacobian;
}

// 从普通单元到标准单元的雅克比矩阵
double TmpEle::Global_to_Local_jacobian(const Point& lp, const std::vector<Point>& gv) const
{
	// lp只是提供维度
	double jacobian = AREA(Pnt[0],Pnt[1],Pnt[2]) / AREA(gv[0],gv[1],gv[2]);
	return jacobian;
}

// 从标准单元到普通单元的雅克比矩阵
std::vector<double> TmpEle::Local_to_Global_jacobian(const std::vector<Point>& lp, const std::vector<Point>& gv) const
{
	std::vector<double> gj(lp.size());
	double larea = AREA(Pnt[0], Pnt[1], Pnt[2]);
	double garea = AREA(gv[0], gv[1], gv[2]);
	for (int i = 0; i < gj.size(); i++){
		gj[i] = garea / larea;
	}
	return gj;
}

// 从普通单元到标准单元的雅克比矩阵
std::vector<double> TmpEle::Global_to_Local_jacobian(const std::vector<Point>& lp, const std::vector<Point>& gv) const
{
	std::vector<double> gj(lp.size());
	double larea = AREA(Pnt[0], Pnt[1], Pnt[2]);
	double garea = AREA(gv[0], gv[1], gv[2]);
	for (int i = 0; i < gj.size(); i++){
		gj[i] = larea / garea;
	}
	return gj;
}
<file_sep>/Heat方程-C++显格式/TmpEle.h
/*************************************************************************
	> File Name: TmpEle.h
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时06分26秒
 ************************************************************************/

#ifndef TMPELE_H
#define TMPELE_H

#include "Point.h"
class TmpEle
{
    friend std::ostream& operator<< (std::ostream&, const TmpEle&);

private:
    std::vector<Point> Pnt;
    std::vector<Point> GaussPnt;
	std::vector<double> GaussWeight;

public:
	// 构造函数
	TmpEle();
	TmpEle(int);
	void buildTE(int);

	Point getPnt(int);
	double getArea();
	std::vector<Point> getGaussPnt();
	std::vector<double> getGaussWeight();
    // 标准单元到普通单元
	Point Local_to_Global(const Point&, const std::vector<Point> &) const;
	// 普通单元到标准单元
	Point Global_to_Local(const Point&, const std::vector<Point> &) const;
	// 标准单元到普通单元
	std::vector<Point> Local_to_Global(const std::vector<Point> &,const std::vector<Point> & ) const;
	// 普通单元到标准单元
	std::vector<Point> Global_to_Local(const std::vector<Point>&, const std::vector<Point>&) const;
	// 从标准单元到普通单元的雅克比矩阵
	double Local_to_Global_jacobian(const Point&, const std::vector<Point>&) const;
	// 从普通单元到标准单元的雅克比矩阵
	double Global_to_Local_jacobian(const Point&, const std::vector<Point>&) const;
	// 从标准单元到普通单元的雅克比矩阵
	std::vector<double> Local_to_Global_jacobian(const std::vector<Point>&, const std::vector<Point>&) const;
	// 从普通单元到标准单元的雅克比矩阵
	std::vector<double> Global_to_Local_jacobian(const std::vector<Point>&, const std::vector<Point>&) const; 
};

#endif
<file_sep>/Heat方程-C++显格式/Point.cpp
/*************************************************************************
	> File Name: Point.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时01分06秒
 ************************************************************************/
#include "Point.h"

//非成员接口函数
std::ostream& operator<< (std::ostream& os, const Point& p)
{
	for (auto i : p.x)
		os << i << '\t';
	return os;
}

std::istream& operator>> (std::istream& is, Point& p)
{
	for (auto &i : p.x)
		is >> i;
	return is;
}

Point operator+ (const Point& p1, const Point& p2)
{
	Point p;
	for (int i = 0; i != 2; ++i)
		p.x[i] = p1.x[i] + p2.x[i];
	return p;
}

Point operator- (const Point& p1, const Point& p2)
{
	Point p;
	for (int i = 0; i != 2; ++i)
		p.x[i] = p1.x[i] - p2.x[i];
	return p;
}

Point midpoint(const Point& p1, const Point& p2)
{
	Point p;
	for (int i = 0; i != 2; ++i)
		p.x[i] = (p1.x[i] + p2.x[i]) / 2.0;
	return p;
}

double distance(const Point& p1, const Point& p2)
{
	double p = 0.;
	for (int i = 0; i != 2; ++i)
		p += (p1.x[i] - p2.x[i]) * (p1.x[i] - p2.x[i]);
	return sqrt(p);
}

Point barycenter(const std::vector<Point>& p, const double * w)
{
	double bc[2] = {0,0};
	int k = p.size();
	if (w == NULL){
		for (int i = 0; i < k; i++){
			bc[0] += p[i][0];
			bc[1] += p[i][1];
		}
		bc[0] /= k;
		bc[1] /= k;
	}
	else{
		double sw = 0;
		for (int i = 0; i < k; i++) sw += w[i];
		for (int i = 0; i < k; i++){
			bc[0] += w[i] * p[i][0];
			bc[1] += w[i] * p[i][1];
		}
		bc[0] /= sw;
		bc[1] /= sw;
	}
	return bc;
}

double AREA(const Point& p1, const Point& p2, const Point& p3)
{
	double area = ((p2.x[0] - p1.x[0]) * (p3.x[1] - p1.x[1]) - 
		(p3.x[0] - p1.x[0]) * (p2.x[1] - p1.x[1])) * 0.5;
	if (area < 0)
		area = -area;
	return area;
}
// 构造函数
Point::Point()
{
	for (int i = 0; i != 2; ++i)
		x[i] = 0;

}

Point::Point(const double* data)
{
	for (int i = 0; i != 2; ++i)
		x[i] = data[i];
}

Point::Point(const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] = p.x[i];
}

Point::Point(const double& v1, const double& v2)
{
	x[0] = v1;
	x[1] = v2;
}
// 成员函数
Point& Point::operator= (const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] = p.x[i];
	return *this;
}

double& Point::operator[] (int i)
{
	return x[i];
}

const double& Point::operator[] (int i) const
{
	return x[i];
}

double Point::length() const
{
	double v = 0.0;
	for (auto i : x)
		v += i * i;
	return sqrt(v);
}

Point& Point::operator+= (const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] += p.x[i];
	return *this;
}

Point& Point::operator-= (const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] -= p.x[i];
	return *this;
}

Point& Point::operator*= (const double& s)
{
	for (auto &i : x)
		i *= s;
	return *this;
}

Point& Point::operator/= (const double& s)
{
	for (auto &i : x)
		i /= s;
	return *this;
}
<file_sep>/README.md
# Finite element algorithm
Matlab;C++
<file_sep>/Heat方程-C++隐格式-非零边界/PDE.h
/*************************************************************************
	> File Name: PDE.h
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时42分05秒
 ************************************************************************/

#ifndef PDE_H
#define PDE_H

#include "Point.h"

const double Pi = 4.0 * atan(1.0);

class PDE
{	
public:
	// 偏微分方程的边界条件(第一类边界条件)
	double u_boundary(const Point &, double);
	// 偏微分方程的右端项
	double f(const Point &, double);
	// 偏微分方程解析解
	double u_exact(const Point&, double);

    // 偏微分方程解析解的梯度
    std::vector<double> u_exact_grad(const Point &, double);
};

#endif
<file_sep>/Possion方程--C++(CVT自适应)/Matrix.cpp
#include "Matrix.h"

#define PI (4.0*atan(1.0))

////////////////////////////////////////////////////////////////////////////////
// 类Point
////////////////////////////////////////////////////////////////////////////////

//非成员接口函数
std::ostream& operator<< (std::ostream& os, const Point& p)
{
	for (auto i : p.x)
		os << i << '\t';
	return os;
}

std::istream& operator>> (std::istream& is, Point& p)
{
	for (auto &i : p.x)
		is >> i;
	return is;
}

Point operator+ (const Point& p1, const Point& p2)
{
	Point p;
	for (int i = 0; i != 2; ++i)
		p.x[i] = p1.x[i] + p2.x[i];
	return p;
}

Point operator- (const Point& p1, const Point& p2)
{
	Point p;
	for (int i = 0; i != 2; ++i)
		p.x[i] = p1.x[i] - p2.x[i];
	return p;
}

Point midpoint(const Point& p1, const Point& p2)
{
	Point p;
	for (int i = 0; i != 2; ++i)
		p.x[i] = (p1.x[i] + p2.x[i]) / 2.0;
	return p;
}

double distance(const Point& p1, const Point& p2)
{
	double p = 0.;
	for (int i = 0; i != 2; ++i)
		p += (p1.x[i] - p2.x[i]) * (p1.x[i] - p2.x[i]);
	return sqrt(p);
}

Point barycenter(const std::vector<Point>& p, const double * w)
{
	double bc[2] = {0,0};
	int k = p.size();
	if (w == NULL){
		for (int i = 0; i < k; i++){
			bc[0] += p[i][0];
			bc[1] += p[i][1];
		}
		bc[0] /= k;
		bc[1] /= k;
	}
	else{
		double sw = 0;
		for (int i = 0; i < k; i++) sw += w[i];
		for (int i = 0; i < k; i++){
			bc[0] += w[i] * p[i][0];
			bc[1] += w[i] * p[i][1];
		}
		bc[0] /= sw;
		bc[1] /= sw;
	}
	return bc;
}

double AREA(const Point& p1, const Point& p2, const Point& p3)
{
	double area = ((p2.x[0] - p1.x[0]) * (p3.x[1] - p1.x[1]) - 
		(p3.x[0] - p1.x[0]) * (p2.x[1] - p1.x[1])) * 0.5;
	if (area < 0)
		area = -area;
	return area;
}
// 构造函数
Point::Point()
{
	for (int i = 0; i != 2; ++i)
		x[i] = 0;

}

Point::Point(const double* data)
{
	for (int i = 0; i != 2; ++i)
		x[i] = data[i];
}

Point::Point(const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] = p.x[i];
}

Point::Point(const double& v1, const double& v2)
{
	x[0] = v1;
	x[1] = v2;
}
// 成员函数
Point& Point::operator= (const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] = p.x[i];
	return *this;
}

double& Point::operator[] (int i)
{
	return x[i];
}

const double& Point::operator[] (int i) const
{
	return x[i];
}

double Point::length() const
{
	double v = 0.0;
	for (auto i : x)
		v += i * i;
	return sqrt(v);
}

Point& Point::operator+= (const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] += p.x[i];
	return *this;
}

Point& Point::operator-= (const Point& p)
{
	for (int i = 0; i != 2; ++i)
		x[i] -= p.x[i];
	return *this;
}

Point& Point::operator*= (const double& s)
{
	for (auto &i : x)
		i *= s;
	return *this;
}

Point& Point::operator/= (const double& s)
{
	for (auto &i : x)
		i /= s;
	return *this;
}

////////////////////////////////////////////////////////////////////////////////
//类TmpEle
////////////////////////////////////////////////////////////////////////////////
// 非成员函数
std::ostream& operator<< (std::ostream& os, const TmpEle& te)
{	
	for (auto i : te.Pnt)
		os << i << '\n';
	for (auto i : te.GaussPnt)
		os << i << '\n';
	for (auto i : te.GaussWeight)
		os << i << '\n';
	return os;
}

// 构造函数
TmpEle::TmpEle()
{
	Pnt.resize(3);
	Pnt[0][0]=0.0;
	Pnt[0][1]=0.0; 
	Pnt[1][0]=1.0;
	Pnt[1][1]=0.0; 
	Pnt[2][0]=0.0;
	Pnt[2][1]=1.0; 
	
	GaussPnt.resize(7);
	GaussWeight.resize(7);
	GaussPnt[0][0]=0.333333333333;
	GaussPnt[1][0]=0.470142064105;
	GaussPnt[2][0]=0.470142064105;
	GaussPnt[3][0]=0.059715871790;
	GaussPnt[4][0]=0.101286507323;
	GaussPnt[5][0]=0.101286507323;
	GaussPnt[6][0]=0.797426985353;
	GaussPnt[0][1]=0.333333333333;
	GaussPnt[1][1]=0.470142064105;
	GaussPnt[2][1]=0.059715871790;
	GaussPnt[3][1]=0.470142064105;
	GaussPnt[4][1]=0.101286507323;
	GaussPnt[5][1]=0.797426985353;
	GaussPnt[6][1]=0.101286507323;
	GaussWeight[0]=0.225000000000;
	GaussWeight[1]=0.132394152788;
	GaussWeight[2]=0.132394152789;
	GaussWeight[3]=0.132394152788;
	GaussWeight[4]=0.125939180545;
	GaussWeight[5]=0.125939180545; 
	GaussWeight[6]=0.125939180545; 
}

TmpEle::TmpEle(int i)
{
  Pnt.resize(3);
  Pnt[0][0]=0.0;
  Pnt[0][1]=0.0; 
  Pnt[1][0]=1.0;
  Pnt[1][1]=0.0; 
  Pnt[2][0]=0.0;
  Pnt[2][1]=1.0; 
  switch(i){
     case 1:
        GaussPnt.resize(1);
		GaussWeight.resize(1);
		GaussPnt[0][0]=0.3333333333333333;
		GaussPnt[0][1]=0.3333333333333333;
		GaussWeight[0]=1.0;
        break;
     case 2:
        GaussPnt.resize(3);
		GaussWeight.resize(3);
		GaussPnt[0][0]=0.166666666667;
		GaussPnt[1][0]=0.166666666667;
		GaussPnt[2][0]=0.666666666667;
		GaussPnt[0][1]=0.166666666667;
		GaussPnt[1][1]=0.666666666667;
		GaussPnt[2][1]=0.166666666667;
		GaussWeight[0]=0.333333333333;
		GaussWeight[1]=0.333333333334;
		GaussWeight[2]=0.333333333333;
        break;
     case 3:
        GaussPnt.resize(4);
		GaussWeight.resize(4);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.200000000000;
		GaussPnt[2][0]=0.200000000000;
		GaussPnt[3][0]=0.600000000000;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.200000000000;
		GaussPnt[2][1]=0.600000000000;
		GaussPnt[3][1]=0.200000000000;
		GaussWeight[0]=-0.562500000000;
		GaussWeight[1]=0.520833333333;
		GaussWeight[2]=0.520833333334;
		GaussWeight[3]=0.520833333333;
        break;
     case 4:
        GaussPnt.resize(6);
		GaussWeight.resize(6);
		GaussPnt[0][0]=0.445948490916;
		GaussPnt[1][0]=0.445948490916;
		GaussPnt[2][0]=0.108103018168;
		GaussPnt[3][0]=0.091576213510;
		GaussPnt[4][0]=0.091576213510;
		GaussPnt[5][0]=0.816847572980;
		GaussPnt[0][1]=0.445948490916;
		GaussPnt[1][1]=0.108103018168;
		GaussPnt[2][1]=0.445948490916;
		GaussPnt[3][1]=0.091576213510;
		GaussPnt[4][1]=0.816847572980;
		GaussPnt[5][1]=0.091576213510;
		GaussWeight[0]=0.223381589678;
		GaussWeight[1]=0.223381589678;
		GaussWeight[2]=0.223381589678;
		GaussWeight[3]=0.109951743655;
		GaussWeight[4]=0.109951743656;
		GaussWeight[5]=0.109951743655;
        break;
     case 5:
        GaussPnt.resize(7);
		GaussWeight.resize(7);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.470142064105;
		GaussPnt[2][0]=0.470142064105;
		GaussPnt[3][0]=0.059715871790;
		GaussPnt[4][0]=0.101286507323;
		GaussPnt[5][0]=0.101286507323;
		GaussPnt[6][0]=0.797426985353;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.470142064105;
		GaussPnt[2][1]=0.059715871790;
		GaussPnt[3][1]=0.470142064105;
		GaussPnt[4][1]=0.101286507323;
		GaussPnt[5][1]=0.797426985353;
		GaussPnt[6][1]=0.101286507323;
		GaussWeight[0]=0.225000000000;
		GaussWeight[1]=0.132394152788;
		GaussWeight[2]=0.132394152789;
		GaussWeight[3]=0.132394152788;
		GaussWeight[4]=0.125939180545;
		GaussWeight[5]=0.125939180545; 
		GaussWeight[6]=0.125939180545; 
        break;
     case 6:
        GaussPnt.resize(12);
		GaussWeight.resize(12);
		GaussPnt[0][0]=0.249286745171;
		GaussPnt[1][0]=0.249286745171;
		GaussPnt[2][0]=0.501426509658;
		GaussPnt[3][0]=0.063089014492;
		GaussPnt[4][0]=0.063089014492;
		GaussPnt[5][0]=0.873821971017;
		GaussPnt[6][0]=0.310352451034;
		GaussPnt[7][0]=0.636502499121;
		GaussPnt[8][0]=0.053145049845;
		GaussPnt[9][0]=0.636502499121;
		GaussPnt[10][0]=0.310352451034;
		GaussPnt[11][0]=0.053145049845;
		GaussPnt[0][1]=0.249286745171;
		GaussPnt[1][1]=0.501426509658;
		GaussPnt[2][1]=0.249286745171;
		GaussPnt[3][1]=0.063089014492;
		GaussPnt[4][1]=0.873821971017;
		GaussPnt[5][1]=0.063089014492;
		GaussPnt[6][1]=0.636502499121;
		GaussPnt[7][1]=0.053145049845;
		GaussPnt[8][1]=0.310352451034;
		GaussPnt[9][1]=0.310352451034;
		GaussPnt[10][1]=0.053145049845;
		GaussPnt[11][1]=0.636502499121;
		GaussWeight[0]=0.116786275726;
		GaussWeight[1]=0.116786275726;
		GaussWeight[2]=0.116786275726;
		GaussWeight[3]=0.050844906370;
		GaussWeight[4]=0.050844906370;
		GaussWeight[5]=0.050844906370; 
		GaussWeight[6]=0.082851075618;
		GaussWeight[7]=0.082851075618; 
		GaussWeight[8]=0.082851075618; 
		GaussWeight[9]=0.082851075618; 
		GaussWeight[10]=0.082851075618; 
		GaussWeight[11]=0.082851075618; 
        break;
     case 7:
        GaussPnt.resize(13);
		GaussWeight.resize(13);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.260345966079;
		GaussPnt[2][0]=0.260345966079;
		GaussPnt[3][0]=0.479308067842;
		GaussPnt[4][0]=0.065130102902;
		GaussPnt[5][0]=0.065130102902;
		GaussPnt[6][0]=0.869739794196;
		GaussPnt[7][0]=0.312865496005;
		GaussPnt[8][0]=0.638444188570;
		GaussPnt[9][0]=0.048690315425;
		GaussPnt[10][0]=0.638444188570;
		GaussPnt[11][0]=0.312865496005;
		GaussPnt[12][0]=0.048690315425;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.260345966079;
		GaussPnt[2][1]=0.479308067842;
		GaussPnt[3][1]=0.260345966079;
		GaussPnt[4][1]=0.065130102902;
		GaussPnt[5][1]=0.869739794196;
		GaussPnt[6][1]=0.065130102902;
		GaussPnt[7][1]=0.638444188570;
		GaussPnt[8][1]=0.048690315425;
		GaussPnt[9][1]=0.312865496005;
		GaussPnt[10][1]=0.312865496005;
		GaussPnt[11][1]=0.048690315425;
		GaussPnt[12][1]=0.638444188570;
		GaussWeight[0]=-0.149570044468;
		GaussWeight[1]=0.175615257433;
		GaussWeight[2]=0.175615257433;
		GaussWeight[3]=0.175615257433;
		GaussWeight[4]=0.053347235609;
		GaussWeight[5]=0.053347235609; 
		GaussWeight[6]=0.053347235609;
		GaussWeight[7]=0.077113760890; 
		GaussWeight[8]=0.077113760890; 
		GaussWeight[9]=0.077113760890; 
		GaussWeight[10]=0.077113760890; 
		GaussWeight[11]=0.077113760890; 
		GaussWeight[12]=0.077113760890;  
        break;
     case 10:
        GaussPnt.resize(25);
		GaussWeight.resize(25);
		GaussPnt[0][0]=0.333333333333;
		GaussPnt[1][0]=0.485577633384;
		GaussPnt[2][0]=0.485577633384;
		GaussPnt[3][0]=0.028844733233;
		GaussPnt[4][0]=0.109481575485;
		GaussPnt[5][0]=0.109481575485;
		GaussPnt[6][0]=0.781036849030;
		GaussPnt[7][0]=0.307939838764;
		GaussPnt[8][0]=0.550352941821;
		GaussPnt[9][0]=0.141707219415;
		GaussPnt[10][0]=0.550352941821;
		GaussPnt[11][0]=0.307939838764;
		GaussPnt[12][0]=0.141707219415;
		GaussPnt[13][0]=0.246672560640;
		GaussPnt[14][0]=0.728323904597;
		GaussPnt[15][0]=0.025003534763;
		GaussPnt[16][0]=0.728323904597;
		GaussPnt[17][0]=0.246672560640;
		GaussPnt[18][0]=0.025003534763;
		GaussPnt[19][0]=0.066803251012;
		GaussPnt[20][0]=0.923655933587;
		GaussPnt[21][0]=0.009540815400;
		GaussPnt[22][0]=0.923655933587;
		GaussPnt[23][0]=0.066803251012;
		GaussPnt[24][0]=0.009540815400;
		GaussPnt[0][1]=0.333333333333;
		GaussPnt[1][1]=0.485577633384;
		GaussPnt[2][1]=0.028844733233;
		GaussPnt[3][1]=0.485577633384;
		GaussPnt[4][1]=0.109481575485;
		GaussPnt[5][1]=0.781036849030;
		GaussPnt[6][1]=0.109481575485;
		GaussPnt[7][1]=0.550352941821;
		GaussPnt[8][1]=0.141707219415;
		GaussPnt[9][1]=0.307939838764;
		GaussPnt[10][1]=0.307939838764;
		GaussPnt[11][1]=0.141707219415;
		GaussPnt[12][1]=0.550352941821;
		GaussPnt[13][1]=0.728323904597;
		GaussPnt[14][1]=0.025003534763;
		GaussPnt[15][1]=0.246672560640;
		GaussPnt[16][1]=0.246672560640;
		GaussPnt[17][1]=0.025003534763;
		GaussPnt[18][1]=0.728323904597;
		GaussPnt[19][1]=0.923655933587;
		GaussPnt[20][1]=0.009540815400;
		GaussPnt[21][1]=0.066803251012;
		GaussPnt[22][1]=0.066803251012;
		GaussPnt[23][1]=0.009540815400;
		GaussPnt[24][1]=0.923655933587;
		GaussWeight[0]=0.090817990383;
		GaussWeight[1]=0.036725957756;
		GaussWeight[2]=0.036725957756;
		GaussWeight[3]=0.036725957756;
		GaussWeight[4]=0.045321059436;
		GaussWeight[5]=0.045321059436; 
		GaussWeight[6]=0.045321059436;
		GaussWeight[7]=0.072757916845; 
		GaussWeight[8]=0.072757916845; 
		GaussWeight[9]=0.072757916845; 
		GaussWeight[10]=0.072757916845; 
		GaussWeight[11]=0.072757916845; 
		GaussWeight[12]=0.072757916845; 
		GaussWeight[13]=0.028327242531; 
		GaussWeight[14]=0.028327242531; 
		GaussWeight[15]=0.028327242531; 
		GaussWeight[16]=0.028327242531; 
		GaussWeight[17]=0.028327242531; 
		GaussWeight[18]=0.028327242531; 
		GaussWeight[19]=0.009421666964; 
		GaussWeight[20]=0.009421666964; 
		GaussWeight[21]=0.009421666964; 
		GaussWeight[22]=0.009421666964; 
		GaussWeight[23]=0.009421666964; 
		GaussWeight[24]=0.009421666964; 
        break;
     case 12:
        GaussPnt.resize(33);
		GaussWeight.resize(33);
		GaussPnt[0][0]=0.488217389774;
		GaussPnt[1][0]=0.488217389774;
		GaussPnt[2][0]=0.023565220452;
		GaussPnt[3][0]=0.439724392294;
		GaussPnt[4][0]=0.439724392294;
		GaussPnt[5][0]=0.120551215411;
		GaussPnt[6][0]=0.271210385012;
		GaussPnt[7][0]=0.271210385012;
		GaussPnt[8][0]=0.457579229976;
		GaussPnt[9][0]=0.127576145542;
		GaussPnt[10][0]=0.127576145542;
		GaussPnt[11][0]=0.744847708917;
		GaussPnt[12][0]=0.021317350453;
		GaussPnt[13][0]=0.021317350453;
		GaussPnt[14][0]=0.957365299094;
		GaussPnt[15][0]=0.275713269686;
		GaussPnt[16][0]=0.608943235780;
		GaussPnt[17][0]=0.115343494535;
		GaussPnt[18][0]=0.608943235780;
		GaussPnt[19][0]=0.275713269686;
		GaussPnt[20][0]=0.115343494535;
		GaussPnt[21][0]=0.281325580990;
		GaussPnt[22][0]=0.695836086788;
		GaussPnt[23][0]=0.022838332222;
		GaussPnt[24][0]=0.695836086788;
		GaussPnt[25][0]=0.281325580990;
		GaussPnt[26][0]=0.022838332222;
		GaussPnt[27][0]=0.116251915908;
		GaussPnt[28][0]=0.858014033544;
		GaussPnt[29][0]=0.025734050548;
		GaussPnt[30][0]=0.858014033544;
		GaussPnt[31][0]=0.116251915908;
		GaussPnt[32][0]=0.025734050548;
		GaussPnt[0][1]=0.488217389774;
		GaussPnt[1][1]=0.023565220452;
		GaussPnt[2][1]=0.488217389774;
		GaussPnt[3][1]=0.439724392294;
		GaussPnt[4][1]=0.120551215411;
		GaussPnt[5][1]=0.439724392294;
		GaussPnt[6][1]=0.271210385012;
		GaussPnt[7][1]=0.457579229976;
		GaussPnt[8][1]=0.271210385012;
		GaussPnt[9][1]=0.127576145542;
		GaussPnt[10][1]=0.744847708917;
		GaussPnt[11][1]=0.127576145542;
		GaussPnt[12][1]=0.021317350453;
		GaussPnt[13][1]=0.957365299094;
		GaussPnt[14][1]=0.021317350453;
		GaussPnt[15][1]=0.608943235780;
		GaussPnt[16][1]=0.115343494535;
		GaussPnt[17][1]=0.275713269686;
		GaussPnt[18][1]=0.275713269686;
		GaussPnt[19][1]=0.115343494535;
		GaussPnt[20][1]=0.608943235780;
		GaussPnt[21][1]=0.695836086788;
		GaussPnt[22][1]=0.022838332222;
		GaussPnt[23][1]=0.281325580990;
		GaussPnt[24][1]=0.281325580990;
		GaussPnt[25][1]=0.022838332222;
		GaussPnt[26][1]=0.695836086788;
		GaussPnt[27][1]=0.858014033544;
		GaussPnt[28][1]=0.025734050548;
		GaussPnt[29][1]=0.116251915908;
		GaussPnt[30][1]=0.116251915908;
		GaussPnt[31][1]=0.025734050548;
		GaussPnt[32][1]=0.858014033544;
		GaussWeight[0]=0.025731066440;
		GaussWeight[1]=0.025731066440;
		GaussWeight[2]=0.025731066440;
		GaussWeight[3]=0.043692544538;
		GaussWeight[4]=0.043692544538;
		GaussWeight[5]=0.043692544538; 
		GaussWeight[6]=0.062858224218;
		GaussWeight[7]=0.062858224218; 
		GaussWeight[8]=0.062858224218; 
		GaussWeight[9]=0.034796112931; 
		GaussWeight[10]=0.034796112931; 
		GaussWeight[11]=0.034796112931; 
		GaussWeight[12]=0.006166261052; 
		GaussWeight[13]=0.006166261052; 
		GaussWeight[14]=0.006166261052; 
		GaussWeight[15]=0.040371557766; 
		GaussWeight[16]=0.040371557766; 
		GaussWeight[17]=0.040371557766; 
		GaussWeight[18]=0.040371557766; 
		GaussWeight[19]=0.040371557766; 
		GaussWeight[20]=0.040371557766; 
		GaussWeight[21]=0.022356773202; 
		GaussWeight[22]=0.022356773202; 
		GaussWeight[23]=0.022356773202; 
		GaussWeight[24]=0.022356773202; 
		GaussWeight[25]=0.022356773202; 
		GaussWeight[26]=0.022356773202; 
		GaussWeight[27]=0.017316231109; 
		GaussWeight[28]=0.017316231109; 
		GaussWeight[29]=0.017316231109; 
		GaussWeight[30]=0.017316231109; 
		GaussWeight[31]=0.017316231109; 
		GaussWeight[32]=0.017316231109; 
        break;
     default:
        std::cout<<"There are no Gauss quadrature points with accuracy of "<<i<<"-th order"<<std::endl;
  }
}

// 成员函数

void TmpEle::buildTE(int i)
{
	TmpEle tmpEle(i);
	for (auto x : tmpEle.Pnt)
		Pnt.push_back(x);
	for (auto x : tmpEle.GaussPnt)
		GaussPnt.push_back(x);
	for (auto x : tmpEle.GaussWeight)
		GaussWeight.push_back(x);
}

Point TmpEle::getPnt(int i)
{
	return Pnt[i];
}

std::vector<Point> TmpEle::getGaussPnt()
{
	return GaussPnt;
}

std::vector<double> TmpEle::getGaussWeight()
{
	return GaussWeight;
}

double TmpEle::getArea()
{
	return 0.5*((Pnt[1][0] - Pnt[0][0]) * (Pnt[2][1] - Pnt[0][1])
				- (Pnt[1][1] - Pnt[0][1]) * (Pnt[2][0] - Pnt[0][0]));
}

// 标准单元到普通单元
Point TmpEle::Local_to_Global(const Point & lp, const std::vector<Point> & gv) const
{
	Point gp;
	double lambda[3];
	double area = AREA(Pnt[0], Pnt[1], Pnt[2]);
	lambda[0] = AREA(lp, Pnt[1], Pnt[2]) / area;
	lambda[1] = AREA(lp, Pnt[2], Pnt[0]) / area;
	lambda[2] = AREA(lp, Pnt[0], Pnt[1]) / area;
	gp[0] = lambda[0] * gv[0][0] + lambda[1] * gv[1][0] + lambda[2] * gv[2][0];
	gp[1] = lambda[0] * gv[0][1] + lambda[1] * gv[1][1] + lambda[2] * gv[2][1];
	return gp;
}

// 普通单元到标准单元
Point TmpEle::Global_to_Local(const Point& gp, const std::vector<Point>& gv) const
{
	Point lp;
	double area = AREA(gv[0],gv[1],gv[2]);

	lp[0] = ((gv[2][1] - gv[0][1]) * (gp[0] - gv[0][0]) - 
		(gv[2][0] - gv[0][0]) * (gp[1] - gv[0][1])) / (2 * area);
	lp[1] = (-(gv[1][1] - gv[0][1]) * (gp[0] - gv[0][0]) +
		(gv[1][0] - gv[0][0]) * (gp[1] - gv[0][1])) / (2 * area);

	return lp;
}

// 标准单元到普通单元
std::vector<Point> TmpEle::Local_to_Global(const std::vector<Point>& lp, const std::vector<Point> & gv) const
{
	std::vector<Point> gp(lp.size());
	double area = AREA(Pnt[0], Pnt[1], Pnt[2]);
	for(int i=0; i<gp.size(); i++){
		double lambda[3];
		lambda[0] = AREA(lp[i], Pnt[1], Pnt[2]) / area;
		lambda[1] = AREA(lp[i], Pnt[2], Pnt[0]) / area;
		lambda[2] = AREA(lp[i], Pnt[0], Pnt[1]) / area;
		gp[i][0] = lambda[0] * gv[0][0] + lambda[1] * gv[1][0] + lambda[2] * gv[2][0];
		gp[i][1] = lambda[0] * gv[0][1] + lambda[1] * gv[1][1] + lambda[2] * gv[2][1];
	}
	return gp;
}

// 普通单元到标准单元
std::vector<Point> TmpEle::Global_to_Local(const std::vector<Point>& gp, 
	const std::vector<Point>& gv) const
{
	std::vector<Point> lp;
	double area = AREA(gv[0],gv[1],gv[2]);

	for (auto point : gp){
		Point p;
		p[0] = ((gv[2][1] - gv[0][1]) * (point[0] - gv[0][0]) - 
				(gv[2][0] - gv[0][0]) * (point[1] - gv[0][1])) / (2 * area);
		p[1] = (-(gv[1][1] - gv[0][1]) * (point[0] - gv[0][0]) + 
				(gv[1][0] - gv[0][0]) * (point[1] - gv[0][1])) / (2 * area);
		lp.push_back(p);
	}
	return lp;
}

// 从标准单元到普通单元的雅克比矩阵
double TmpEle::Local_to_Global_jacobian(const Point& lp, const std::vector<Point>& gv) const
{
	// lp只是提供维度
	double jacobian =  AREA(gv[0],gv[1],gv[2]) / AREA(Pnt[0], Pnt[1], Pnt[2]);
	return jacobian;
}

// 从普通单元到标准单元的雅克比矩阵
double TmpEle::Global_to_Local_jacobian(const Point& lp, const std::vector<Point>& gv) const
{
	// lp只是提供维度
	double jacobian = AREA(Pnt[0],Pnt[1],Pnt[2]) / AREA(gv[0],gv[1],gv[2]);
	return jacobian;
}

// 从标准单元到普通单元的雅克比矩阵
std::vector<double> TmpEle::Local_to_Global_jacobian(const std::vector<Point>& lp, const std::vector<Point>& gv) const
{
	std::vector<double> gj(lp.size());
	double larea = AREA(Pnt[0], Pnt[1], Pnt[2]);
	double garea = AREA(gv[0], gv[1], gv[2]);
	for (int i = 0; i < gj.size(); i++){
		gj[i] = garea / larea;
	}
	return gj;
}

// 从普通单元到标准单元的雅克比矩阵
std::vector<double> TmpEle::Global_to_Local_jacobian(const std::vector<Point>& lp, const std::vector<Point>& gv) const
{
	std::vector<double> gj(lp.size());
	double larea = AREA(Pnt[0], Pnt[1], Pnt[2]);
	double garea = AREA(gv[0], gv[1], gv[2]);
	for (int i = 0; i < gj.size(); i++){
		gj[i] = larea / garea;
	}
	return gj;
}


////////////////////////////////////////////////////////////////////////////////
//类Mesh
////////////////////////////////////////////////////////////////////////////////

// 非成员函数
std::ostream& operator<< (std::ostream& os, const Mesh& M)
{
	for (auto i : M.Pnt)
		os << i << '\n';

	for (auto &i : M.Ele){
		for (auto j : i)
			os << j << '\t';
		os << '\n';
	}

	for (auto i : M.BndPnt)
		os << i << '\n';
		os << '\n';
	return os;
}

// 成员函数
void Mesh::readData(const std::string& f)
{
	int i;
	std::ifstream is(f,std::ifstream::in); // 以只读的模式打开
	is >> i;
	Pnt.resize(i);
	for (int j = 0; j < i; j++)
		is >> Pnt[j][0] >> Pnt[j][1];

	int n;
	is >> n;
	Ele.resize(n);
	for (int j = 0; j < n; j++){
		Ele[j].resize(3);
		is >> Ele[j][0] >> Ele[j][1] >> Ele[j][2];
	}

	int bn;
	is >> bn;
	BndPnt.resize(bn);
	for (int j = 0; j < bn; j++)
		is >> BndPnt[j];

	is.close();
}

int Mesh::getEleVtx(int i, int j)
{
	return Ele[i][j];
}

std::vector<int> Mesh::getEleVtx(int i)
{
	return Ele[i];
}

Point Mesh::getPnt(int i)
{
	return Pnt[i];
}

std::vector<Point> Mesh::getPnt(std::vector<int>& vt)
{
	std::vector<Point> vec;
	for (int x : vt){
		Point point = Pnt[x];
		vec.push_back(point);
	}
	return vec;
}

std::vector<int> Mesh::getBndPnt()
{
	return BndPnt;
}

int Mesh::n_element()
{
	return Ele.size();
}

int Mesh::n_point()
{
	return Pnt.size();
}

int Mesh::n_boundaryPoint()
{
	return BndPnt.size();
}

void Mesh::neighbor(std::vector<std::vector<int>>& neighbor)
{
    auto n_ele = n_element();   // 得到剖分小三角形的数量
    // 初始化
    neighbor.resize(n_ele);
    for (decltype(neighbor.size()) i = 0; i != n_ele; ++i){
        neighbor[i].resize(3);
    }
    for (decltype(neighbor.size()) i = 0; i != n_ele; ++i){
        // 检测顶点所对应的边是不是边界边
        bool boundary_mark_00 = false, boundary_mark_01 = false;
        bool boundary_mark_10 = false, boundary_mark_11 = false;
        bool boundary_mark_20 = false, boundary_mark_21 = false;
        for (auto boundary : BndPnt){
            // 第一个顶点所对应的边是不是边界边
            if (boundary == Ele[i][1]) boundary_mark_00 = true;
            if (boundary == Ele[i][2]) boundary_mark_01 = true;
            // 第二个顶点所对应的边是不是边界边
            if (boundary == Ele[i][0]) boundary_mark_10 = true;
            if (boundary == Ele[i][2]) boundary_mark_11 = true;
            // 第三个顶点所对应的边是不是边界边
            if (boundary == Ele[i][0]) boundary_mark_20 = true;
            if (boundary == Ele[i][1]) boundary_mark_21 = true;
        }
        // 如果是边界边的话邻居单元就是自己
        if (boundary_mark_00 && boundary_mark_01) neighbor[i][0] = i;
        if (boundary_mark_10 && boundary_mark_11) neighbor[i][1] = i;
        if (boundary_mark_20 && boundary_mark_21) neighbor[i][2] = i;
        
        for (decltype(neighbor.size()) j = 0; j != n_ele; ++j){
            if (i == j) continue;
            bool mark_00 = false, mark_01 = false;
            bool mark_10 = false, mark_11 = false;
            bool mark_20 = false, mark_21 = false;
            for (decltype(neighbor.size()) k = 0; k != 3; ++k){
                // 第一个顶点对应边的邻居单元(如果是边界边的话邻居单元就是自己)
                if (Ele[j][k] == Ele[i][1]) mark_00 = true;
                if (Ele[j][k] == Ele[i][2]) mark_01 = true;
                // 第二个顶点对应边的邻居单元(如果是边界边的话邻居单元就是自己)
                if (Ele[j][k] == Ele[i][0]) mark_10 = true;
                if (Ele[j][k] == Ele[i][2]) mark_11 = true;
                // 第三个顶点对应边的邻居单元(如果是边界边的话邻居单元就是自己)
                if (Ele[j][k] == Ele[i][0]) mark_20 = true;
                if (Ele[j][k] == Ele[i][1]) mark_21 = true;
            }
            if (mark_00 && mark_01) neighbor[i][0] = j;
            if (mark_10 && mark_11) neighbor[i][1] = j;
            if (mark_20 && mark_21) neighbor[i][2] = j; 
        }
        
    }
    
}


////////////////////////////////////////////////////////////////////////////////
//偏微分方程,边界条件和初值条件结构体,PDE
//
////////////////////////////////////////////////////////////////////////////////

// 边界条件(第一类边界条件)

double PDE::u_boundary(const Point& p)
{
	return 0;
}

double PDE::f(const Point& p)
{
	double value;
	value = - 2 * PI * PI * exp(PI * (p[0] + p[1])) * (sin(PI * p[0]) * cos(PI * p[1]) + cos(PI * p[0]) * sin(PI * p[1]));
	return value; 
}

double PDE::u_exact(const Point& p)
{
	double value;
	value = exp(PI * (p[0] + p[1])) * sin(PI * p[0]) * sin(PI * p[1]);
	return value;
}

////////////////////////////////////////////////////////////////////////////////
//类Matrix
//
////////////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////////////


void Matrix::GaussSeidel(const std::vector<std::vector<double>>& A_matrix, std::vector<double>& u, 
       const std::vector<double>& rhs)
{
   std::vector<double> last_u(A_matrix.size());  
   do{
      last_u=u; 
      double error=0.0;     
      for(int i=0;i<A_matrix.size();i++){
         double temp=rhs[i];
         for(int j=0; j<A_matrix[i].size(); j++){
             if(j != i){
                temp-=A_matrix[i][j]*u[j];
             }   
         }
         u[i]=temp/A_matrix[i][i];
      }
      u_int Vsize=u.size();
      for(int k=0;k<Vsize;k++){
         error+=(last_u[k]-u[k])*(last_u[k]-u[k]);
      }
      error=sqrt(error);
      
      if(error<1.0e-10)
          break;
   }while(1); 
}

void Matrix::SolveLinearASparse(const Eigen::SparseMatrix<double>& A_matrix_sparse,std::vector<double>& u,
       const std::vector<double>& rhs, int k)
{
    auto b_size = rhs.size();
    Eigen::VectorXd b;
    b.resize(b_size);
    for (size_t i = 0; i != b_size; ++i)
        b(i) = rhs[i];
    Eigen::VectorXd u_h_now;
    // 求解器声明
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IdentityPreconditioner> bicg;
    Eigen::GMRES<Eigen::SparseMatrix<double>, Eigen::IdentityPreconditioner> gmres;
    Eigen::DGMRES<Eigen::SparseMatrix<double>,Eigen::IdentityPreconditioner> dgmres;
    Eigen::MINRES<Eigen::SparseMatrix<double>,Eigen::Lower|Eigen::Upper,Eigen::IdentityPreconditioner> minres;
    switch (k){
        // 迭代法
        case 1:
            cg.compute(A_matrix_sparse);
            u_h_now = cg.solve(b);
            std::cout << "CG:\t #iterations: " << cg.iterations() << ", estimated error: " << cg.error() 
                    << std::endl;
            break;
        case 2:
            bicg.compute(A_matrix_sparse);
            u_h_now = bicg.solve(b);
            std::cout << "BiCGSTAB:\t #iterations: " << bicg.iterations() << ", estimated error: " 
                    << bicg.error() << std::endl;
            break;
        case 3:
            gmres.compute(A_matrix_sparse);
            u_h_now = gmres.solve(b);
            std::cout << "GMRES:\t #iterator: " << gmres.iterations() << ", estimated error: "
                    << gmres.error() << std::endl;
            break;
        case 4:
            dgmres.compute(A_matrix_sparse);
            u_h_now = dgmres.solve(b);
            std::cout << "DGMRES:\t #iterator: " << dgmres.iterations() << ", estimate error: "
                    << dgmres.error() << std::endl;
            break;
        case 5:
            minres.compute(A_matrix_sparse);
            u_h_now = minres.solve(b);
            std::cout << "MINRES:\t #iterator: " << minres.iterations() << ", estimate error: "
                    << minres.error() << std::endl;
            break;
    }

    // 将求解结果存储到u中
    for (size_t i = 0; i != u.size(); ++i)
        u[i] = u_h_now(i);
}

Matrix::Matrix(const std::string& file, int i)
{
	mesh.readData(file);
	tmpEle.buildTE(i);
}

// 成员函数
// 基函数的构造(面积坐标)
void Matrix::basisValue(const Point& p, const std::vector<Point>& v, std::vector<double>& val)
{
	val.resize(3); // 分片线性单元
	double area = AREA(v[0], v[1], v[2]);
	val[0] = AREA(p,v[1],v[2]);
	val[1] = AREA(v[0],p,v[2]);
	val[2] = AREA(v[0],v[1],p);
	val[0] /= area;
	val[1] /= area;
	val[2] /= area;
}

void Matrix::basisValue(const std::vector<Point>& p, const std::vector<Point>& v, 
	std::vector<std::vector<double>>& val)
{
	int n = p.size();
	val.resize(3);
	for(int i = 0; i < 3; i++)
		val[i].resize(n);
	double area = AREA(v[0],v[1],v[2]);

	for (int i=0; i < n; i++){
		val[0][i] = AREA(p[i],v[1],v[2]);
		val[1][i] = AREA(v[0],p[i],v[2]);
		val[2][i] = AREA(v[0],v[1],p[i]);
		val[0][i] /= area;
		val[1][i] /= area;
		val[2][i] /= area;
	}
}

void Matrix::basisGrad(const Point&, const std::vector<Point>& v, 
	std::vector<std::vector<double>>& val)
{
	val.resize(3);
	val[0].resize(2);
	val[1].resize(2);
	val[2].resize(2);

	double area = AREA(v[0],v[1],v[2]) * 2;  // 面积坐标到直角坐标的雅克比为2s

	val[0][0] = (v[1][1] - v[2][1]) / area;
	val[0][1] = (v[2][0] - v[1][0]) / area;

	val[1][0] = (v[2][1] - v[0][1]) / area;
	val[1][1] = (v[0][0] - v[2][0]) / area;

	val[2][0] = (v[0][1] - v[1][1]) / area;
	val[2][1] = (v[1][0] - v[0][0]) / area;
}

void Matrix::basisGrad(const std::vector<Point>& p, const std::vector<Point>& v, 
	std::vector<std::vector<std::vector<double>>>& val)
{
	int n = p.size();
	val.resize(3);
	for (int i = 0; i < 3; i++){
		val[i].resize(n);
	}
	for (int i = 0; i < 3; i++){
		for (int j = 0; j < n; j++){
			val[i][j].resize(2);
		}
	}
	double area = AREA(v[0],v[1],v[2]) * 2; 
	for (int i = 0; i < n; i++){

		val[0][i][0] = (v[1][1] - v[2][1]) / area;
		val[0][i][1] = (v[2][0] - v[1][0]) / area;

		val[1][i][0] = (v[2][1] - v[0][1]) / area;
		val[1][i][1] = (v[0][0] - v[2][0]) / area;

		val[2][i][0] = (v[0][1] - v[1][1]) / area;
		val[2][i][1] = (v[1][0] - v[0][0]) / area;
	}
}

void Matrix::buildEleMatrix_RHS(int n, std::vector<std::vector<double>> &Ele_matrix, std::vector<double> &Ele_Rhs)
{
	Ele_matrix.resize(3);  // 单元刚度矩阵(左端),维度是3*3的
	Ele_Rhs.resize(3);   // 单元刚度矩阵(右端),维度是3*1的
	for(int i=0; i<3; i++){
		Ele_matrix[i].resize(3);
	}

	// 得到第n个单元的边
	std::vector<int> EV=mesh.getEleVtx(n);
	// 得到单元的三个点
	std::vector<Point> EP(3);
	EP[0]=mesh.getPnt(EV[0]);
	EP[1]=mesh.getPnt(EV[1]);
	EP[2]=mesh.getPnt(EV[2]);
	// 标准三角单元的面积
	double volume = tmpEle.getArea();
	// 得到标准三角单元的Gauss点
	std::vector<Point> GP=tmpEle.getGaussPnt();
	// 得到标准三角单元的Gauss权重
	std::vector<double> GW=tmpEle.getGaussWeight();
	// 把标准单元的Gauss点投影到普通单元上
	std::vector<Point> q_point = tmpEle.Local_to_Global(GP, EP); 
	// 得到标准单元到普通单元的雅克比矩阵
	std::vector<double> jacobian = tmpEle.Local_to_Global_jacobian(GP, EP);
	// 得到三角单元上基函数的值
	std::vector<std::vector<double>> basis_value;
	basisValue(q_point, EP, basis_value);
	// 得到三角单元上梯度算子的值
	std::vector<std::vector<std::vector<double>>> basis_grad;
	basisGrad(q_point, EP, basis_grad);

	// 计算右端的单元刚度矩阵
	for(int i = 0; i < q_point.size(); i++){
		double Jxw = GW[i] * jacobian[i] * volume;
		double f_value = pde.f(q_point[i]);
		for(int j = 0; j < basis_value.size(); j++){
			Ele_Rhs[j] += Jxw*f_value*basis_value[j][i];
		}
	}

	for(int i = 0; i < q_point.size(); i++){
		double Jxw = GW[i] * jacobian[i] * volume;
		for(int k=0; k<basis_grad.size(); k++){
			for(int j=0; j<basis_grad.size(); j++){
				Ele_matrix[j][k] += Jxw * (basis_grad[j][i][0] * basis_grad[k][i][0] + basis_grad[j][i][1] *
                               basis_grad[k][i][1]);
			}
		}
	}
}

void Matrix::buildMatrix_RHS_Sparse(Eigen::SparseMatrix<double>& A_matrix_sparse,std::vector<double>& Rhs)
{
    auto n_ele = mesh.n_element();
    auto n_pnt = mesh.n_point();
    A_matrix_sparse.resize(n_pnt,n_pnt);
    Rhs.resize(n_pnt);
    // 定义一个用于存储稀疏矩阵的三元组
    std::vector<Eigen::Triplet<double>> triple;
    for (size_t i = 0; i != n_ele; ++i){
    	std::vector<int> NV = mesh.getEleVtx(i);

        std::vector<std::vector<double>> Ele_matrix(3);
        for (size_t j = 0; j != 3; ++j){
            Ele_matrix[j].resize(3);
        }
        std::vector<double> Ele_Rhs(3);
        buildEleMatrix_RHS(i,Ele_matrix,Ele_Rhs);
        // 稀疏矩阵存储技术
        // 利用经典的三元组插入方式来存储A_matrix
        triple.push_back(Eigen::Triplet<double>(NV[0],NV[0],Ele_matrix[0][0]));
        triple.push_back(Eigen::Triplet<double>(NV[0],NV[1],Ele_matrix[0][1]));
        triple.push_back(Eigen::Triplet<double>(NV[0],NV[2],Ele_matrix[0][2]));
        triple.push_back(Eigen::Triplet<double>(NV[1],NV[0],Ele_matrix[1][0]));
        triple.push_back(Eigen::Triplet<double>(NV[1],NV[1],Ele_matrix[1][1]));
        triple.push_back(Eigen::Triplet<double>(NV[1],NV[2],Ele_matrix[1][2]));
        triple.push_back(Eigen::Triplet<double>(NV[2],NV[0],Ele_matrix[2][0]));
        triple.push_back(Eigen::Triplet<double>(NV[2],NV[1],Ele_matrix[2][1]));
        triple.push_back(Eigen::Triplet<double>(NV[2],NV[2],Ele_matrix[2][2]));

        // 不利用稀疏矩阵存储右端项Rhs
        Rhs[NV[0]] += Ele_Rhs[0];
        Rhs[NV[1]] += Ele_Rhs[1];
        Rhs[NV[2]] += Ele_Rhs[2];
    }
    // 把三元组转换成稀疏矩阵
    A_matrix_sparse.setFromTriplets(triple.begin(),triple.end());
}

// 对于稀疏矩阵A_matrix_sparse的边界条件处理(dirichlet边界条件)
void Matrix::DealBoundary_Sparse(Eigen::SparseMatrix<double>& A_matrix_sparse, const std::vector<double>& uh,
        std::vector<double>& Rhs)
{
    auto n_pnt = mesh.n_point();
    auto n_bp = mesh.n_boundaryPoint();
    std::vector<int> BV = mesh.getBndPnt(); // 得到边界点
    for (size_t i = 0; i != n_bp; ++i){
        Point Pnt_i = mesh.getPnt(BV[i]); // 得到第i个单元的点
        double val = pde.u_boundary(Pnt_i);
        Rhs[BV[i]] = A_matrix_sparse.coeffRef(BV[i],BV[i]) * val;
        for (size_t j = 0; j != n_pnt; ++j){
            if (j != BV[i]){
                A_matrix_sparse.coeffRef(BV[i],j) = 0.0;
            }
        }
        for (size_t j = 0; j != n_pnt; ++j){
            if (j != BV[i]){
                Rhs[j] -= A_matrix_sparse.coeffRef(j,BV[i]) * val;
                A_matrix_sparse.coeffRef(j,BV[i]) = 0.0;
            }
        }
    }
}

double Matrix::ComputerL2Error(const std::vector<double> &f)
{
	double err=0.0;
	int n_ele=mesh.n_element();
	for(int j=0; j<n_ele; j++){
		std::vector<int> NV = mesh.getEleVtx(j);
		std::vector<Point> EP(3);
		EP[0]=mesh.getPnt(NV[0]);
		EP[1]=mesh.getPnt(NV[1]);
		EP[2]=mesh.getPnt(NV[2]);
		double volume = tmpEle.getArea();
		std::vector<Point> GP = tmpEle.getGaussPnt();

		std::vector<double> GW = tmpEle.getGaussWeight();
		std::vector<Point> q_point = tmpEle.Local_to_Global(GP, EP); 
		std::vector<double> jacobian = tmpEle.Local_to_Global_jacobian(GP, EP);
		std::vector<std::vector<double> > basis_value;
		basisValue(q_point, EP, basis_value);

		for(int i = 0; i < q_point.size(); i++){
			double Jxw = GW[i]*jacobian[i]*volume;
			double f_value = pde.u_exact(q_point[i]);
			double u_h_val = f[NV[0]]*basis_value[0][i]+f[NV[1]]*basis_value[1][i]+f[NV[2]]*basis_value[2][i];
			double df_val = f_value-u_h_val;
			err += Jxw*df_val*df_val;
		}
	}
	err=sqrt(err);
	return err;
}

void Matrix::GradBasis(std::vector<std::vector<std::vector<double>>>& gradbasis)
{
    auto n_ele = mesh.n_element();  // 三角单元的数量
    auto n_point =  mesh.n_point();    // 三角单元点的个数
    gradbasis.resize(3);
    for (decltype(gradbasis.size()) i = 0; i != 3; ++i)
        gradbasis[i].resize(n_ele);
    for (decltype(gradbasis.size()) i = 0; i != 3; ++i)
        for (decltype(gradbasis.size()) j = 0; j != n_ele; ++j)
            gradbasis[i][j].resize(2);
    
    std::vector<int> EV;    // 存储三角单元的边
    std::vector<Point> EP(3);   // 存储三角单元的三个顶点
    for (decltype(gradbasis.size()) i = 0; i != n_ele; ++i){
        // 得到第i个单元的边
        EV = mesh.getEleVtx(i);
        // 得到单元的三个顶点
        EP[0] = mesh.getPnt(EV[0]);
        EP[1] = mesh.getPnt(EV[1]);
        EP[2] = mesh.getPnt(EV[2]);
       
        gradbasis[0][i][0] = (EP[1][1] - EP[2][1]) / (2 * AREA(EP[0],EP[1],EP[2]));
        gradbasis[0][i][1] = (EP[2][0] - EP[1][0]) / (2 * AREA(EP[0],EP[1],EP[2]));

        gradbasis[1][i][0] = (EP[2][1] - EP[0][1]) / (2 * AREA(EP[0],EP[1],EP[2]));
        gradbasis[1][i][1] = (EP[0][0] - EP[2][0]) / (2 * AREA(EP[0],EP[1],EP[2]));

        gradbasis[2][i][0] = (EP[0][1] - EP[1][1]) / (2 * AREA(EP[0],EP[1],EP[2]));
        gradbasis[2][i][1] = (EP[1][0] - EP[0][0]) / (2 * AREA(EP[0],EP[1],EP[2]));
    }
    
}

void Matrix::Gradu(const std::vector<std::vector<std::vector<double>>>& gradbasis, 
        const std::vector<double>& u_h, std::vector<std::vector<double>>& Du)
{
    auto n_ele = mesh.n_element();  // 三角形单元的数量
    Du.resize(n_ele);
    for (decltype(Du.size()) i = 0; i != n_ele; ++i)
        Du[i].resize(2);

    std::vector<int> EV;    // 存储三角单元的边
    std::vector<Point> EP(3);   // 存储三角单元的三个顶点
    
    for (decltype(Du.size()) i =0; i != n_ele; ++i){
        // 得到第i个单元的边
        EV = mesh.getEleVtx(i);
        // 得到单元的三个顶点
        EP[0] = mesh.getPnt(EV[0]);
        EP[1] = mesh.getPnt(EV[1]);
        EP[2] = mesh.getPnt(EV[2]);

        Du[i][0] = u_h[EV[0]] * gradbasis[0][i][0] + u_h[EV[1]] * gradbasis[1][i][0] + u_h[EV[2]] * gradbasis[2][i][0];
        Du[i][1] = u_h[EV[0]] * gradbasis[0][i][1] + u_h[EV[1]] * gradbasis[1][i][1] + u_h[EV[2]] * gradbasis[2][i][1];
    }
    
}

void Matrix::eleVector(std::vector<std::vector<std::vector<double>>>& ve)
{
    // 初始化
    auto n_ele = mesh.n_element();
    ve.resize(n_ele);
    for (decltype(ve.size()) i = 0; i != n_ele; ++i){
        ve[i].resize(2);
        for (decltype(ve.size()) j = 0; j != 2; ++j){
            ve[i][j].resize(3);
        }
    }
    // 边的向量
    for (decltype(ve.size()) i = 0; i != n_ele; ++i){
        for (decltype(ve.size()) j = 0; j != 2; ++j){
            // 得到第i个单元的边
            std::vector<int> EV = mesh.getEleVtx(i);
            std::vector<Point> EP(3);
            EP[0] = mesh.getPnt(EV[0]);
            EP[1] = mesh.getPnt(EV[1]);
            EP[2] = mesh.getPnt(EV[2]);
            ve[i][j][0] = EP[2][j] - EP[1][j];
            ve[i][j][1] = EP[0][j] - EP[2][j];
            ve[i][j][2] = EP[1][j] - EP[0][j];
        }
    }
}

void Matrix::normalVector(const std::vector<std::vector<std::vector<double>>>& ve, 
        std::vector<std::vector<std::vector<double>>>& ne)
{
    // 初始化
    auto n_ele = mesh.n_element();
    ne.resize(n_ele);
    for (decltype(ne.size()) i = 0; i != n_ele; ++i){
        ne[i].resize(2);
        for (decltype(ne.size()) j = 0; j != 2; ++j){
            ne[i][j].resize(3);
        }
    }
    // 把ve旋转90度得到ne
    for (decltype(ne.size()) i = 0; i != n_ele; ++i){
        ne[i][0][0] = ve[i][1][0]; ne[i][1][0] = -ve[i][0][0];
        ne[i][0][1] = ve[i][1][1]; ne[i][1][1] = -ve[i][0][1];
        ne[i][0][2] = ve[i][1][2]; ne[i][1][2] = -ve[i][0][2];
    }
}

double Matrix::diameter(const std::vector<Point>& EP)
{
    // 求出三角单元三条边的边长
    auto distance_0 = distance(EP[0],EP[1]);
    auto distance_1 = distance(EP[1],EP[2]);
    auto distance_2 = distance(EP[0],EP[2]);

    double max = distance_0;
    if (max < distance_1)
        max = distance_1;
    if (max < distance_2);
        max = distance_2;
    return max;
}


void Matrix::Estimateresidual(std::vector<double>& u_h,std::vector<double>& eta)
{
    auto n_ele = mesh.n_element();  // 三角形单元的数量
    auto n_point = mesh.n_point();  // 三角形点的个数
    // 得到剖分上每一个单元的梯度
    std::vector<std::vector<std::vector<double>>> gradbasis;
    GradBasis(gradbasis);
    // 得到u_h的梯度
    std::vector<std::vector<double>> Du;
    Gradu(gradbasis,u_h,Du);
    // 得到邻居单元的索引
    std::vector<std::vector<int>> neighbor_index;
    mesh.neighbor(neighbor_index);
    // 得到剖分单元上每一条边的向量
    std::vector<std::vector<std::vector<double>>> ve;     
    eleVector(ve);
    // 得到向量的法向量
    std::vector<std::vector<std::vector<double>>> ne;
    normalVector(ve,ne);
    // 对于边界边的跳量为0
    std::vector<double> edgeJump;
    edgeJump.resize(n_ele);
    for (decltype(edgeJump.size()) i = 0; i != n_ele; ++i){
        auto index_0 = neighbor_index[i][0];
        auto index_1 = neighbor_index[i][1];
        auto index_2 = neighbor_index[i][2];

        std::vector<int> EV = mesh.getEleVtx(i);
        std::vector<Point> EP(3);
        EP[0] = mesh.getPnt(EV[0]);
        EP[1] = mesh.getPnt(EV[1]);
        EP[2] = mesh.getPnt(EV[2]);

        edgeJump[i] = pow(((Du[i][0] - Du[index_0][0]) * ne[i][0][0] + (Du[i][1] - Du[index_0][1]) * ne[i][1][0]),2.0) 
            +  pow(((Du[i][0] - Du[index_1][0]) * ne[i][0][1] + (Du[i][1] - Du[index_1][1]) * ne[i][1][1]),2.0)  + 
            pow(((Du[i][0] - Du[index_2][0]) * ne[i][0][2] + (Du[i][1] - Du[index_2][1]) * ne[i][1][2]),2.0);
    }
    // 初始化
    eta.resize(n_ele);
    // 计算残差||r||_{L^2}^2
    for (decltype(eta.size()) i = 0; i != n_ele; ++i ){
        std::vector<int> NV = mesh.getEleVtx(i);
        std::vector<Point> EP(3);
        EP[0] = mesh.getPnt(NV[0]);
        EP[1] = mesh.getPnt(NV[1]);
        EP[2] = mesh.getPnt(NV[2]);
        double volume = tmpEle.getArea();
        // 得到高斯点
        std::vector<Point> GP = tmpEle.getGaussPnt();
        // 得到高斯权重
        std::vector<double> GW = tmpEle.getGaussWeight();
        // 把高斯点投影到普通单元
        std::vector<Point> q_point = tmpEle.Local_to_Global(GP,EP);
        std::vector<double> jacobian = tmpEle.Local_to_Global_jacobian(GP,EP);
        // 得到改三角单元的直径
        double dm = diameter(EP);
        for (decltype(eta.size()) j = 0; j != q_point.size(); ++j){
            double Jxw = GW[j] * jacobian[j] * volume;
            double f_value = pde.f(q_point[j]);
            eta[i] += Jxw * f_value * f_value * dm * dm;
        }
        // 加入跳量
        eta[i] += 0.5 * dm * edgeJump[i];
    }
}









<file_sep>/Possion方程--C++(CVT自适应)/density.cpp
/*************************************************************************
	> File Name: density.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年08月20日 星期一 22时33分51秒
 ************************************************************************/

#include "density.h"

void density(const std::vector<Point>& node, const std::vector<std::vector<int>>& ele, const std::vector<double>& eta)
{
    auto N = node.size();   // 网格节点数目
    auto NT = ele.size();   // 单元数目
    // 生成文件"dens_nodes.dat"
    std::ofstream out("dens_nodes.dat",std::ofstream::out);     // 以只读的模式打开

    auto radi = 6371000;
    std::vector<std::vector<double>> sp_crd;
    // 初始化
    sp_crd.resize(N);
    for (decltype(sp_crd.size()) i = 0; i != N; ++i)
        sp_crd[i].resize(3);
    
    for (decltype(sp_crd.size()) i = 0; i != N; ++i){
        auto x = 0.5 * node[i][0] / radi;
        auto y = 0.5 * node[i][1] / radi;
        auto xy2 = x * x + y * y;
        sp_crd[i][0] = radi * 2 * x / (1+xy2);
        sp_crd[i][1] = radi * 2 * y / (1+xy2);
        sp_crd[i][2] = radi * (1-xy2) / (1+xy2);
    }
    // 输出到文档中
    for (decltype(sp_crd.size()) i = 0; i != sp_crd.size(); ++i){
        out << sp_crd[i][0] << "\t" << sp_crd[i][1] << "\t" << sp_crd[i][2] << "\n";
    }
    out.close();
    
    // 生成文件"dens_valus.dat"
    // 检查节点i是否在单元j内,是为1,不是为0
    std::vector<std::vector<int>> p2t;
    p2t.resize(N);
    for (decltype(p2t.size()) i = 0; i != N; ++i)
        p2t[i].resize(NT);

    for (decltype(p2t.size()) i = 0; i != N; ++i){
        for (decltype(p2t.size()) j = 0; j != NT; ++j){
            bool direction_0 = false;
            bool direction_1 = false;
            bool direction_2 = false;
            if (i == ele[j][0]) direction_0 = true;
            if (i == ele[j][1]) direction_1 = true;
            if (i == ele[j][2]) direction_2 = true;
            if (direction_0 || direction_1 || direction_2)
                p2t[i][j] = 1;
            else
                p2t[i][j] = 0;
        }
    }
    
    // 生成单元尺寸
    std::vector<double> rho;    // 点的尺寸N
    std::vector<double> area;
    area.resize(NT);    // 初始化
    rho.resize(N);    // 初始化

    std::vector<double> temporary_rho;
    temporary_rho.resize(NT);   // 边的尺寸

    for (size_t i = 0; i != NT; ++i){
        Point point_0 = node[ele[i][0]];
        Point point_1 = node[ele[i][1]];
        Point point_2 = node[ele[i][2]];
        area[i] = AREA(point_0,point_1,point_2);
        auto distance_0 = distance(point_0,point_1);
        auto distance_1 = distance(point_0,point_2);
        auto distance_2 = distance(point_1,point_2);
        auto htri = (distance_0 + distance_1 + distance_2) / 3;
        // 计算暂时的rho
        temporary_rho[i] = pow(eta[i],2.0) / pow(htri,4.0);
    }
    // 转换为点的(这里有一个尺寸的变化)
    for (decltype(p2t.size()) i = 0; i != N; ++i){
        double temporary_0 = 0.0;
        double temporary_1 = 0.0;
        for (decltype(p2t[i].size()) j = 0; j != NT; ++j){
            temporary_0 += p2t[i][j] * area[j] * temporary_rho[j];
            temporary_1 += p2t[i][j] * area[j];
        }
        rho[i] = temporary_0 / temporary_1;
    }

    for (size_t i = 0; i != 5; ++i){
        for (size_t j = 0; j != NT; ++j){
            double temporary = 0.0;
            for (size_t k = 0; k != N; ++k){
                temporary += p2t[k][j] * rho[k] / 3.0;
            }
            temporary_rho[j] = temporary;
        }
        for (size_t j = 0; j != N; ++j){
            double temporary_0 = 0.0;
            double temporary_1 = 0.0;
            for (size_t k = 0; k != NT; ++k){
                temporary_0 += p2t[j][k] * temporary_rho[k] * area[k];
                temporary_1 += p2t[j][k] * area[k];
            }
            rho[j] = temporary_0 / temporary_1; 
        }
    }

    // 写入文件
    out.open("dens_valus.dat",std::ofstream::out);
    for (decltype(rho.size()) i = 0; i != rho.size(); ++i){
        out << rho[i] << "\n";
    }
    out.close();
}


<file_sep>/Heat方程-C++隐格式/Makefile
c++ = g++

VERSION = -std=c++0x

all: Test

Test: Test.o Point.o Mesh.o Parabolic.o TmpEle.o PDE.o
	$(c++) $(VERSION) -o Test Test.o Point.o Mesh.o Parabolic.o TmpEle.o PDE.o

Test.o: Test.cpp
	$(c++) $(VERSION) -c Test.cpp

Point.o: Point.cpp Point.h
	$(c++) $(VERSION) -c Point.cpp

Mesh.o: Mesh.cpp Mesh.h
	$(c++) $(VERSION) -c Mesh.cpp

Parabolic.o: Parabolic.cpp Parabolic.h
	$(c++) $(VERSION) -c Parabolic.cpp

TmpEle.o: TmpEle.cpp TmpEle.h
	$(c++) $(VERSION) -c TmpEle.cpp

PDE.o: PDE.cpp PDE.h
	$(c++) $(VERSION) -c PDE.cpp
<file_sep>/Possion方程--C++/Test.cpp
#include "Matrix.h"

using std::cout;
using std::endl;
using std::cin;

void PointTest()
{
	////////////////////////////////////////////////////////////////////////////
	// 测试Point类
	////////////////////////////////////////////////////////////////////////////
	// 第一种初始化的方式
	cout << "*********************************************************" << endl;
	cout << "\n\t\t类Point的测试\n\n\n";
	cout << "*********************************************************" << endl;
	Point p;
	cout << "默认初始化:\n" << p << endl;

	//第二种初始化的方式
	double a[] = {1,1};
	Point p1(a);
	cout << "初始化方式p(数组):\n" << p1 << endl;

	//第三种初始化的方式
	Point p2(p1);
	cout << "初始化方式p(Point类型):\n" << p2 << endl;

	//第四种初值化方式
	Point p3 = p2;
	cout << "初始化方式p = Point类型:\n" << p3 << endl;

	Point pp1(0.5,0.5);
	cout << "初始化方式p(double, double):\n" << pp1 << endl;

	cout << "重载运算符[]:\n" << p3[0] << '\t' << p3[1] << endl;

	cout << "输出Point类型的长度:\n " << p3.length() << endl;

	p3 += p2;
	cout << "重载运算符+=:\n" << p3 << endl;

	p3 -= p2;
	cout << "重载运算符-=:\n" << p3 << endl;

	p3 *= 2;
	cout << "重载运算符*=:\n" << p3 << endl;

	p3 /= 2;
	cout << "重载运算符/=:\n" << p3 << endl;

	p3 = p2 + p1;
	cout << "重载运算符+:\n " << p3 << endl;

	p3 = p2 - p1;
	cout << "重载运算符-:\n" << p3 << endl;
	cout << "输入两个数:" << endl;
	cin >> p3;
	cout << "重载输入流 >>:\n" <<p3 << endl;

	cout << "输出两个点的中点:\n" << midpoint(p1,p3) << endl;

	cout << "输出两个点的距离:\n" << distance(p1,p3) << endl;

	Point pp2(0,0), pp3(0.5,0), pp4(0,0.5);
	cout << "输出三角形的面积:\n" << AREA(pp2,pp3,pp4) << endl;

	std::vector<Point> pp{pp2,pp3,pp4};
	Point barycenter0 = barycenter(pp, NULL);
	cout << "输出三角形的质心:\n" << barycenter0 << endl;
}

void TmpEleTest()
{
	////////////////////////////////////////////////////////////////////////////
	// 测试TmpEle类
	////////////////////////////////////////////////////////////////////////////
	cout << "*********************************************************" << endl;
	cout << "\n\t\t类TmpEle的测试\n\n\n";
	cout << "*********************************************************" << endl;
	TmpEle te;
	cout << "默认初始化以及输出流重载的测试:\n" << te <<endl;

	TmpEle te1(1);
	cout << "给私有成员赋予一阶的Gauss点和权重:\n" << te1 << endl;

	TmpEle te2(2);
	cout << "给私有成员赋予二阶的Gauss点和权重:\n" << te2 << endl;

	TmpEle te3(3);
	cout << "给私有成员赋予三阶的Gauss点和权重:\n" << te3 << endl;

	TmpEle te4(4);
	cout << "给私有成员赋予四阶的Gauss点和权重:\n" << te4 << endl;

	TmpEle te5(5);
	cout << "给私有成员赋予五阶的Gauss点和权重:\n" << te5 << endl;

	TmpEle te6(6);
	cout << "给私有成员赋予六阶的Gauss点和权重:\n" << te6 << endl;

	TmpEle te7(7);
	cout << "给私有成员赋予七阶的Gauss点和权重:\n" << te7 << endl;

	TmpEle te10(10);
	cout << "给私有成员赋予十阶的Gauss点和权重:\n" << te10 << endl;

	TmpEle te12(12);
	cout << "给私有成员赋予十二阶的Gauss点和权重：\n" << te12 << endl;


	Point p4 = te12.getPnt(1);
	cout << "成员函数getPnt(int):\n" << p4 << endl;

	double Area = te12.getArea();
	cout << "成员函数getArea():\n" << Area << endl;

	cout << "成员函数getGaussPnt():" << endl;
	std::vector<Point> GP = te12.getGaussPnt();
	for (auto i : GP)
		cout << i << endl;

	cout << "成员函数getGaussWeight():" << endl; 
	std::vector<double> GW = te12.getGaussWeight();
	for (auto i : GW)
		cout << i << endl;

	Point p5(0,0), p6(0.5,0), p7(0,0.5);
	std::vector<Point> v1{p5,p6,p7};

	Point ppp1(0,0);
	Point p8 = te.Local_to_Global(ppp1,v1);
	cout << "测试函数Local_to_Global:\n" << "(0,0)->>\t"<< p8 << endl;
	
	Point ppp2(0,1);
	p8 = te.Local_to_Global(ppp2,v1);
	cout << "(0,1)->>\t" << p8 << endl;

	Point ppp3(1,0);
	p8 = te.Local_to_Global(ppp3,v1);
	cout << "(1,0)->>\t" << p8 << endl;

	Point ppp4(0,0);
	Point p9 = te.Global_to_Local(ppp4,v1);
	cout << "测试函数Global_to_Local:\n" << "(0,0)->> \t" << p9 << endl;

	Point ppp5(0,0.5);
	p9 = te.Global_to_Local(ppp5,v1);
	cout << "(0,0.5)->>\t" << p9 << endl;

	Point ppp6(0.5,0);
	p9 = te.Global_to_Local(ppp6,v1);
	cout << "(0.5,0)->>\t" << p9 << endl;

	std::vector<Point> v2{ppp1,ppp2,ppp3};
	std::vector<Point> vv = te.Local_to_Global(v2,v1);
	cout << "测试函数Local_to_Global:\n";
	for (int i = 0; i < v2.size(); i++)
		cout << v2[i] << "-->>" <<vv[i] << endl;

	std::vector<Point> vv1 = te.Global_to_Local(v1,v1);
	cout << "测试函数Global_to_Local:\n";
	for (int i = 0; i < v1.size(); i++)
		cout << v1[i] << "-->>" << vv1[i] << endl;

	double jacobian = te.Local_to_Global_jacobian(ppp1,v1);
	cout << "从标准单元到普通单元的雅克比矩阵:\n" << jacobian << endl;

	jacobian = te.Global_to_Local_jacobian(ppp1,v1);
	cout << "从普通单元到标准单元的雅克比矩阵:\n" << jacobian << endl;

	std::vector<double> jacobian1 = te.Local_to_Global_jacobian(v1,v1);
	cout << "从标准单元到普通单元的雅克比矩阵:\n";
	for (auto i : jacobian1)
		cout << i << "\t";
	cout << endl;

	jacobian1 = te.Global_to_Local_jacobian(v1,v1);
	cout << "从普通单元到标准单元的雅克比矩阵:\n";
	for (auto i : jacobian1)
		cout << i << "\t";
	cout << endl;
}

void MeshTest()
{
	 ////////////////////////////////////////////////////////////////////////////
    //测试类Mesh
    ////////////////////////////////////////////////////////////////////////////
	cout << "*********************************************************" << endl;
	cout << "\n\t\t类Mesh的测试\n\n\n";
	cout << "*********************************************************" << endl;

	Mesh mesh;
	cout << "测试默认构造函数:\n" << mesh << endl;

	mesh.readData("Trimesh.txt");
	cout << "测试函数readData\n" << mesh << endl;

	int ele = mesh.getEleVtx(2,2);  // 选取第三行边的第三个边
	cout << "测试函数getEleData\n" << ele << endl;

	std::vector<int> ele1 = mesh.getEleVtx(2);
	cout << "测试函数getEleData:\n";
	for (auto i : ele1)
		cout << i << "\t";
	cout << endl;

	Point p = mesh.getPnt(2);
	cout << "测试函数getPnt:\n" << p << endl;

	int n_ele = mesh.n_element();
	cout << "得到三角单元的个数:\n" << n_ele << endl;

	int n_point = mesh.n_point();
	cout << "得到剖分单元点的个数:\n" << n_point << endl;

	int n_boundaryPoint = mesh.n_boundaryPoint();
	cout << "得到剖分单元边界点的个数:\n" << n_boundaryPoint << endl;

	std::vector<int> num{1,2,3,4,5,6,7,8,9};
	std::vector<Point> point = mesh.getPnt(num);
	cout << "测试函数getPnt,输入一系列整数:\n";
	for (auto i : point)
		cout << i << "\n";
	cout << endl;
}

void PDETest()
{
	cout << "*********************************************************" << endl;
	cout << "\n\t\t类PDE的测试\n\n\n";
	cout << "*********************************************************" << endl;

	PDE pde;
	double p[2] = {1,1};
	cout << "PDE中成员函数u_boundary的测试:\n" << pde.u_boundary(p) << endl;
	Point p1(1,1);
	cout << pde.u_boundary(p) << endl;

	cout << "PDE中成员函数u_exact的测试:\n" << pde.u_exact(p) << endl;
	cout << pde.u_exact(p1) << endl;

	cout << "PDE中成员函数f的测试:\n" << pde.f(p) << endl;
	cout << pde.f(p1) << endl;
}


void MatrixTest()
{
	////////////////////////////////////////////////////////////////////////////
	//测试类Possion
	////////////////////////////////////////////////////////////////////////////
	cout << "*********************************************************" << endl;
	cout << "\n\t\t类Matrix的测试\n\n\n";
	cout << "*********************************************************" << endl;

	Matrix possion;
	cout << "****基函数构造测试****" << endl;
	Point p0(0,0);
	Point p1(1,0);
	Point p2(0,1);
	std::vector<Point> Te{p0,p1,p2};
	Point center = barycenter(Te,NULL);
	cout << "输出三角单元的质心:\n" << center << endl;
	std::vector<double> basis;
	cout << "标准单元上基函数的值:\n" << endl;
	possion.basisValue(p0,Te,basis);
	for (auto i : basis)
		cout << i << "\t";
	cout << endl;

	possion.basisValue(p1,Te,basis);
	for (auto i : basis)
		cout << i << "\t";
	cout << endl;

	possion.basisValue(p2,Te,basis);
	for (auto i : basis)
		cout << i << "\t";
	cout << endl;

	cout << "测试标准单元上的基函数的值,(重载basisValue函数):\n";
	std::vector<std::vector<double>> basis1;
	std::vector<Point> Te1{p0,p1,p2,p0,p1,p2};
	possion.basisValue(Te1,Te,basis1);
	for (auto& i : basis1){
		for (auto j : i)
			cout << j << "\t";
		cout << "\n";
	}
	cout << endl;

	cout << "****梯度算子计算测试****" << endl;
	std::vector<std::vector<double>> basisGrad1;
	possion.basisGrad(p0,Te,basisGrad1);
	for (auto &i : basisGrad1){
		for (auto j : i)
			cout << j << "\t";
		cout << "\n";
	}
	cout << endl;

	cout << "重载梯度算子的函数测试:\n";
	std::vector<std::vector<std::vector<double>>> basisGrad2;
	possion.basisGrad(Te,Te,basisGrad2);
	for (auto& i : basisGrad2){
		for (auto& j : i){
			for (auto k : j){
				cout << k << "\t";
			}
			cout << "\n";
		}
	} 
	cout <<endl;

	////////////////////////////////////////////////////////////////////////////
	//左端刚度矩阵和右端矩阵的测试(未经过边界处理)
	////////////////////////////////////////////////////////////////////////////

	Matrix Possion_matrix("Trimesh.txt",7);
	PDE pde;

	int n_pnt = Possion_matrix.mesh.n_point();
	Possion_matrix.u_h.resize(n_pnt);
	
	std::vector<std::vector<double>> A_matrix;
	std::vector<double> Rhs;

	Possion_matrix.buildMatrix_RHS(A_matrix, Rhs);

	cout << "输出(未经过边界处理)左端刚度矩阵和右端项:\n" << endl;
	// 输出左端刚度矩阵和右端项
	for (auto &i : A_matrix){
		for (auto j : i)
			cout << j << "\t";
		cout << "\n";
	}
	cout << endl;

	cout << "\n\n\n";
	for (auto i : Rhs)
		cout << i << "\n";
	cout << endl;

	Possion_matrix.DealBoundary(A_matrix, Possion_matrix.u_h, Rhs);

	cout << "输出(经过边界处理)左端刚度矩阵和右端项:\n" << endl;
	// 输出左端刚度矩阵和右端项
	for (auto &i : A_matrix){
		for (auto j : i)
			cout << j << "\t";
		cout << "\n";
	}
	cout << endl;

	cout << "\n\n\n";
	for (auto i : Rhs)
		cout << i << "\n";
	cout << endl;

 	Possion_matrix.GaussSeidel(A_matrix,Possion_matrix.u_h, Rhs);
	cout << "数值解:\n" << endl;

	for (auto i : Possion_matrix.u_h)
		cout << i << "\n";
	cout << endl;


 	double L2Error=Possion_matrix.ComputerL2Error(Possion_matrix.u_h);
   	std::cout<<"L2Error="<<L2Error<<std::endl;


}

int main()
{
	PointTest();

	TmpEleTest();
	
	MeshTest();	

	PDETest();

	MatrixTest();

	return 0;
}<file_sep>/Heat方程-C++显格式/Parabolic.h
#ifndef PARABOLIC_H
#define PARABOLIC_H

#include "Point.h"
#include "PDE.h"
#include "TmpEle.h"
#include "Mesh.h"

extern double Dt;           // 时间剖分的步长

class Parabolic
{
private:
	// 数据结构
	PDE pde;	
	Mesh mesh;     
	TmpEle tmpEle;  
	std::vector<double> u_h;
    // 私有函数
	// gauss-seidel迭代
	void GaussSeidel(std::vector<std::vector<double>>&,std::vector<double>&, std::vector<double>&);   
    // 利用Eigen库求解稀疏代数方程组
    void SolveLinearASparse(Eigen::SparseMatrix<double>&,std::vector<double>&, std::vector<double>&, int);
	// 基函数的构造
	void basisValue(Point&, std::vector<Point>&, std::vector<double>&);
	void basisValue(std::vector<Point>&, std::vector<Point>&, std::vector<std::vector<double>>&);
	// 梯度算子的计算
	void basisGrad(Point&, std::vector<Point>&, std::vector<std::vector<double>>&);
	void basisGrad(std::vector<Point>&, std::vector<Point>&, 
                    std::vector<std::vector<std::vector<double>>>&);
	// 单元刚度矩阵
	void buildEleMatrix_RHS(int, std::vector<std::vector<double>>&, 
                            std::vector<std::vector<double>>&, std::vector<double>&, double);    
    // 总刚度矩阵(稀疏存储)
    void buildMatrix_RHS_Sparse(Eigen::SparseMatrix<double>&, 
                                Eigen::SparseMatrix<double>&, std::vector<double>&, double); 
    // 稀疏矩阵的边界条件的处理
    void DealBoundary_Sparse(Eigen::SparseMatrix<double>&,
                             std::vector<double>&, double);
    // 求解抛物方程的全离散显格式
    void fullDiscreteFormat(Eigen::SparseMatrix<double> &,
                            const Eigen::SparseMatrix<double>&,
                            const std::vector<double> &, std::vector<double> &);
    // 求解抛物方程的全离散隐格式
    void fullDiscreteImplicitFormat(Eigen::SparseMatrix<double> &,
                                    Eigen::SparseMatrix<double> &,
                                    const std::vector<double> &, std::vector<double> &);
	// 计算L2误差
	double ComputerL2Error(const std::vector<double> &, double); 
    // 计算H1误差
    double ComputerH1Error(const std::vector<double> &, double);
    // 输出边界处的值
    void checkBoundary(const std::vector<double> &, double);
    // 输出解析解和数值解
    void checkValue(const std::vector<double> &, double);
public:
    Parabolic() = default;
	Parabolic(const std::string&, int);
    // 求解方程组(全离散显格式)
    void runFullDiscreteFormat(const std::string &);
    // 求解方程组(全离散隐格式)
    void runFullDiscreteImplicitFormat(const std::string &);
};

#endif
<file_sep>/Possion方程--C++(CVT自适应)/CVT_Adaptive_finite.cpp
/*************************************************************************
	> File Name: CVT_Adaptive_finite.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年08月19日 星期日 22时24分34秒
 ************************************************************************/

#include "Possion.h"
#include "density.h"
int main()
{
    // 参数设置
    double maxN = 2e3, theta = 0.5;
    int maxIt = 50;
    // 初始化网格
    system("./mesh_init -f square.poly -q 1 -a 0.01");
    // 把初始网格信息写入文档中
    std::vector<Point> node;
    std::vector<std::vector<int>> ele;
    std::vector<int> boundary_node;
    writeTriangle(node,ele,boundary_node,"Trimesh.txt");
    // 定义记录L2误差
    std::vector<double> L2Errors;
    // 定义记录点的个数
    std::vector<int> numbers_node;
    numbers_node.push_back(node.size());
    // 开始自适应的过程
    for (int index = 0; index != maxIt; ++index){
        // 求解Possion方程并且得到残量型误差指示子
        std::vector<double> eta;    // 存放误差指示子
        double L2Error = Possion("Trimesh.txt","Results.txt",eta,7);
        L2Errors.push_back(L2Error);
        if (node.size() > maxN) break;  // 当网格点的数量大于maxN时候退出循环
        // 优化和细化网格
        density(node,ele,eta);
        system("./mesh_refi -b 1 -p 0.2");
        system("./mesh_opti -b 1");
        std::vector<int> boundary_node1;
        writeTriangle(node,ele,boundary_node1,"Trimesh.txt");
        numbers_node.push_back(node.size());
    }
    // 将L2范误差输出到文档"L2Error.txt"中
    std::ofstream out("L2Error.txt",std::ofstream::out);   
    for (decltype(L2Errors.size()) i = 0; i != L2Errors.size(); ++i)
        out << L2Errors[i] << "\t" << numbers_node[i] << "\n";
    out.close();
    return 0;
}

<file_sep>/Possion方程--C++(CVT自适应)/Possion.cpp
/*************************************************************************
	> File Name: Possion.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年08月19日 星期日 21时23分39秒
 ************************************************************************/
#include "Possion.h"

double Possion(const std::string filename_input,const std::string filename_output, std::vector<double>& eta,int i)
{
    Matrix possion(filename_input,i);
    PDE pde;

    int n_pnt = possion.mesh.n_point();     // 得到网格点的个数
    possion.u_h.resize(n_pnt);  
    
    std::vector<std::vector<double>> A_matrix;
    std::vector<double> Rhs;
    
    // 定义一个用于存储A_matrix的稀疏矩阵A_matrix_sparse
    Eigen::SparseMatrix<double> A_matrix_sparse(A_matrix.size(),A_matrix.size());
    possion.buildMatrix_RHS_Sparse(A_matrix_sparse,Rhs);

    possion.DealBoundary_Sparse(A_matrix_sparse,possion.u_h,Rhs);
    possion.SolveLinearASparse(A_matrix_sparse,possion.u_h,Rhs,3);  // k 可以取1-5;这里取的3:gamres方法
    double L2Error = possion.ComputerL2Error(possion.u_h);

    // 计算eta误差指示子
    possion.Estimateresidual(possion.u_h,eta);
    // 把数据u_h写入结果文档
    std::ofstream out(filename_output);
    for (auto i : possion.u_h)
        out << i << "\n";
    out.close();
    return L2Error;
}

void nodesf2dat(const std::string datafile, std::vector<Point>& node, std::vector<int>& boundary_node)
{
    std::ifstream is(datafile,std::ifstream::in);   // 以只读模式打开
    int numbers; 
    is >> numbers;
    is.close();
    
    std::ifstream read(datafile,std::ifstream::in);     // 以只读模式打开
    std::vector<std::vector<double>> data;
    // 初始化
    data.resize(numbers + 1);
    for (decltype(data.size()) i = 0; i != numbers+1; ++i){
        data[i].resize(4);
    }
    for (decltype(data.size()) i = 0; i != numbers+1; ++i){
        read >> data[i][0] >> data[i][1] >> data[i][2] >> data[i][3];
    }
    read.close();

    // 初始化
    node.resize(numbers);
    for (decltype(node.size()) i = 0; i != numbers; ++i){
        node[i][0] = data[i+1][1];
        node[i][1] = data[i+1][2];
        if (data[i+1][3] != 0)
            boundary_node.push_back(data[i+1][0]-1);
    }
}
 
void trigsf2dat(const std::string datafile, std::vector<std::vector<int>>& ele)
{
    std::ifstream is(datafile,std::ifstream::in);   // 以只读模式打开
    int numbers;
    is >> numbers;
    is.close();

    std::ifstream read(datafile,std::ifstream::in);     // 以只读模式打开
    std::vector<std::vector<double>> infor;
    // 初始化
    infor.resize(numbers + 1);
    for (decltype(infor.size()) i = 0; i != numbers+1; ++i){
        infor[i].resize(4);
    }
    read >> infor[0][0] >> infor[0][1] >> infor[0][2];
    for (decltype(infor.size()) i = 0; i != numbers; ++i){
        read >> infor[i+1][0] >> infor[i+1][1] >> infor[i+1][2] >> infor[i+1][3];
    }
    read.close();

    // 初始化
    ele.resize(numbers);
    for (decltype(ele.size()) i = 0; i != numbers; ++i){
        ele[i].resize(3);
    }
    for (decltype(ele.size()) i = 0; i != numbers; ++i){
        ele[i][0] = infor[i+1][1] - 1;  // triangle 输出边的编号从1开始
        ele[i][1] = infor[i+1][2] - 1;  
        ele[i][2] = infor[i+1][3] - 1;
    }
}

void writeTriangle(std::vector<Point>& node, std::vector<std::vector<int>>& ele, 
                   std::vector<int>& boundary_node, const std::string output_file)
{
    nodesf2dat("nodes.dat",node,boundary_node);
    trigsf2dat("trigs.dat",ele);
    // 输出内容
    std::ofstream out(output_file,std::ofstream::out);
    out << node.size() << "\n";
    for(decltype(node.size()) i = 0; i != node.size(); ++i)
        out << node[i][0] << "\t" << node[i][1] << "\n";
    out << "\n";
   
    out << ele.size() << "\n";
    for (decltype(ele.size()) i = 0; i != ele.size(); ++i)
        out << ele[i][0] << "\t" << ele[i][1] << "\t" << ele[i][2] << "\n";
    out << "\n";

    out << boundary_node.size() << "\n";
    for (decltype(boundary_node.size()) i = 0; i != boundary_node.size(); ++i)
        out << boundary_node[i] << "\n";
    out << "\n";
    out.close();
}

<file_sep>/Heat方程-C++隐格式-非零边界/Mesh.cpp
/*************************************************************************
	> File Name: Mesh.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时39分19秒
 ************************************************************************/

#include "Mesh.h"

// 非成员函数
std::ostream& operator<< (std::ostream& os, const Mesh& M)
{
	for (auto i : M.Pnt)
		os << i << '\n';

	for (auto &i : M.Ele){
		for (auto j : i)
			os << j << '\t';
		os << '\n';
	}

	for (auto i : M.BndPnt)
		os << i << '\n';
		os << '\n';
	return os;
}

// 成员函数
void Mesh::readData(const std::string& f)
{
	int i;
	std::ifstream is(f,std::ifstream::in); // 以只读的模式打开
	is >> i;
	Pnt.resize(i);
	for (int j = 0; j < i; j++)
		is >> Pnt[j][0] >> Pnt[j][1];

	int n;
	is >> n;
	Ele.resize(n);
	for (int j = 0; j < n; j++){
		Ele[j].resize(3);
		is >> Ele[j][0] >> Ele[j][1] >> Ele[j][2];
	}

	int bn;
	is >> bn;
	BndPnt.resize(bn);
	for (int j = 0; j < bn; j++)
		is >> BndPnt[j];

	is.close();
}

int Mesh::getEleVtx(int i, int j)
{
	return Ele[i][j];
}

std::vector<int> Mesh::getEleVtx(int i)
{
	return Ele[i];
}

Point Mesh::getPnt(int i)
{
	return Pnt[i];
}

std::vector<Point> Mesh::getPnt(std::vector<int>& vt)
{
	std::vector<Point> vec;
	for (int x : vt){
		Point point = Pnt[x];
		vec.push_back(point);
	}
	return vec;
}

std::vector<int> Mesh::getBndPnt()
{
	return BndPnt;
}

int Mesh::n_element()
{
	return Ele.size();
}

int Mesh::n_point()
{
	return Pnt.size();
}

int Mesh::n_boundaryPoint()
{
	return BndPnt.size();
}
<file_sep>/Heat方程-C++显格式/Parabolic.cpp
#include "Parabolic.h"

double Dt = 0.01;

void Parabolic::GaussSeidel(std::vector<std::vector<double> > &A_matrix, std::vector<double> &u, std::vector<double> &rhs)
{
   std::vector<double> last_u(A_matrix.size());  
   do{
      last_u=u; 
      double error=0.0;     
      for(int i=0;i<A_matrix.size();i++){
         double temp=rhs[i];
         for(int j=0; j<A_matrix[i].size(); j++){
             if(j != i){
                temp-=A_matrix[i][j]*u[j];
             }   
         }
         u[i]=temp/A_matrix[i][i];
      }
      u_int Vsize=u.size();
      for(int k=0;k<Vsize;k++){
         error+=(last_u[k]-u[k])*(last_u[k]-u[k]);
      }
      error=sqrt(error);
      
      if(error<1.0e-10)
          break;
   }while(1); 
}

void Parabolic::SolveLinearASparse(Eigen::SparseMatrix<double>& A_matrix_sparse,std::vector<double>& u,
        std::vector<double>& rhs, int k)
{
    auto b_size = rhs.size();
    Eigen::VectorXd b;
    b.resize(b_size);
    for (size_t i = 0; i != b_size; ++i)
        b(i) = rhs[i];
    Eigen::VectorXd u_h_now;
    // 求解器声明
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IdentityPreconditioner> bicg;
    Eigen::GMRES<Eigen::SparseMatrix<double>, Eigen::IdentityPreconditioner> gmres;
    Eigen::DGMRES<Eigen::SparseMatrix<double>,Eigen::IdentityPreconditioner> dgmres;
    Eigen::MINRES<Eigen::SparseMatrix<double>,Eigen::Lower|Eigen::Upper,Eigen::IdentityPreconditioner> minres;
    switch (k){
        // 迭代法
        case 1:
            cg.compute(A_matrix_sparse);
            u_h_now = cg.solve(b);
            std::cout << "CG:\t #iterations: " << cg.iterations() << ", estimated error: " << cg.error() 
                    << std::endl;
            break;
        case 2:
            bicg.compute(A_matrix_sparse);
            u_h_now = bicg.solve(b);
            std::cout << "BiCGSTAB:\t #iterations: " << bicg.iterations() << ", estimated error: " 
                    << bicg.error() << std::endl;
            break;
        case 3:
            gmres.compute(A_matrix_sparse);
            u_h_now = gmres.solve(b);
            std::cout << "GMRES:\t #iterator: " << gmres.iterations() << ", estimated error: "
                    << gmres.error() << std::endl;
            break;
        case 4:
            dgmres.compute(A_matrix_sparse);
            u_h_now = dgmres.solve(b);
            std::cout << "DGMRES:\t #iterator: " << dgmres.iterations() << ", estimate error: "
                    << dgmres.error() << std::endl;
            break;
        case 5:
            minres.compute(A_matrix_sparse);
            u_h_now = minres.solve(b);
            std::cout << "MINRES:\t #iterator: " << minres.iterations() << ", estimate error: "
                    << minres.error() << std::endl;
            break;
    }

    // 将求解结果存储到u中
    for (size_t i = 0; i != u.size(); ++i)
        u[i] = u_h_now(i);
}

Parabolic::Parabolic(const std::string& file, int i)
{
	mesh.readData(file);
	tmpEle.buildTE(i);
}

// 成员函数
// 基函数的构造(面积坐标)
void Parabolic::basisValue(Point& p, std::vector<Point>& v, std::vector<double>& val)
{
	val.resize(3); // 分片线性单元
	double area = AREA(v[0], v[1], v[2]);
	val[0] = AREA(p,v[1],v[2]);
	val[1] = AREA(v[0],p,v[2]);
	val[2] = AREA(v[0],v[1],p);
	val[0] /= area;
	val[1] /= area;
	val[2] /= area;
}

void Parabolic::basisValue(std::vector<Point>& p, std::vector<Point>& v, 
	std::vector<std::vector<double>>& val)
{
	int n = p.size();
	val.resize(3);
	for(int i = 0; i < 3; i++)
		val[i].resize(n);
	double area = AREA(v[0],v[1],v[2]);

	for (int i=0; i < n; i++){
		val[0][i] = AREA(p[i],v[1],v[2]);
		val[1][i] = AREA(v[0],p[i],v[2]);
		val[2][i] = AREA(v[0],v[1],p[i]);
		val[0][i] /= area;
		val[1][i] /= area;
		val[2][i] /= area;
	}
}

void Parabolic::basisGrad(Point&, std::vector<Point>& v, 
	std::vector<std::vector<double>>& val)
{
	val.resize(3);
	val[0].resize(2);
	val[1].resize(2);
	val[2].resize(2);

	double area = AREA(v[0],v[1],v[2]) * 2;  // 面积坐标到直角坐标的雅克比为2s

	val[0][0] = (v[1][1] - v[2][1]) / area;
	val[0][1] = (v[2][0] - v[1][0]) / area;

	val[1][0] = (v[2][1] - v[0][1]) / area;
	val[1][1] = (v[0][0] - v[2][0]) / area;

	val[2][0] = (v[0][1] - v[1][1]) / area;
	val[2][1] = (v[1][0] - v[0][0]) / area;
}

void Parabolic::basisGrad(std::vector<Point>& p, std::vector<Point>& v, 
	std::vector<std::vector<std::vector<double>>>& val)
{
	int n = p.size();
	val.resize(3);
	for (int i = 0; i < 3; i++){
		val[i].resize(n);
	}
	for (int i = 0; i < 3; i++){
		for (int j = 0; j < n; j++){
			val[i][j].resize(2);
		}
	}
	double area = AREA(v[0],v[1],v[2]) * 2; 
	for (int i = 0; i < n; i++){

		val[0][i][0] = (v[1][1] - v[2][1]) / area;
		val[0][i][1] = (v[2][0] - v[1][0]) / area;

		val[1][i][0] = (v[2][1] - v[0][1]) / area;
		val[1][i][1] = (v[0][0] - v[2][0]) / area;

		val[2][i][0] = (v[0][1] - v[1][1]) / area;
		val[2][i][1] = (v[1][0] - v[0][0]) / area;
	}
}

void Parabolic::buildEleMatrix_RHS(int n, std::vector<std::vector<double>> &Ele_Amatrix, 
                                   std::vector<std::vector<double>> &Ele_Mmatrix,
                                   std::vector<double> &Ele_Rhs, double t)
{
	Ele_Amatrix.resize(3);  // 单元刚度矩阵(左端a(u,v)),维度是3*3的
    Ele_Mmatrix.resize(3);  // 单元刚度矩阵(左端(u,v)),维度是3*3的
	Ele_Rhs.resize(3, 0.0);   // 单元刚度矩阵(右端),维度是3*1的
	for(int i=0; i<3; i++){
		Ele_Amatrix[i].resize(3, 0.0);
        Ele_Mmatrix[i].resize(3, 0.0);
	}

	// 得到第n个单元的边
	std::vector<int> EV = mesh.getEleVtx(n);
	// 得到单元的三个点
	std::vector<Point> EP(3);
	EP[0] = mesh.getPnt(EV[0]);
	EP[1] = mesh.getPnt(EV[1]);
	EP[2] = mesh.getPnt(EV[2]);
	// 标准三角单元的面积
	double volume = tmpEle.getArea();
	// 得到标准三角单元的Gauss点
	std::vector<Point> GP = tmpEle.getGaussPnt();
	// 得到标准三角单元的Gauss权重
	std::vector<double> GW = tmpEle.getGaussWeight();
	// 把标准单元的Gauss点投影到普通单元上
	std::vector<Point> q_point = tmpEle.Local_to_Global(GP, EP); 
	// 得到标准单元到普通单元的雅克比矩阵
	std::vector<double> jacobian = tmpEle.Local_to_Global_jacobian(GP, EP);
	// 得到三角单元上基函数的值
	std::vector<std::vector<double>> basis_value;
	basisValue(q_point, EP, basis_value);
	// 得到三角单元上梯度算子的值
	std::vector<std::vector<std::vector<double>>> basis_grad;
	basisGrad(q_point, EP, basis_grad);
	// 计算右端的单元刚度矩阵
	for(int i = 0; i < q_point.size(); i++){
		double Jxw = GW[i] * jacobian[i] * volume;
		double f_value = pde.f(q_point[i],t);
		for(int j = 0; j < basis_value.size(); j++){
			Ele_Rhs[j] += Jxw * f_value * basis_value[j][i];
		}
	}
	for(int i = 0; i < q_point.size(); i++){
		double Jxw = GW[i] * jacobian[i] * volume;
		for(int k = 0; k<basis_grad.size(); k++){
			for(int j = 0; j<basis_grad.size(); j++){
				// 左端的单元刚度矩阵(a(u,v))
                Ele_Amatrix[j][k] += Jxw * (basis_grad[j][i][0] * 
                                            basis_grad[k][i][0] + basis_grad[j][i][1] * basis_grad[k][i][1]);
            }
		}
        for (int k = 0; k < basis_value.size(); ++k){
            for (int j = 0; j < basis_value.size(); ++j){
                // 左端的单元刚度矩阵((u,v))
                Ele_Mmatrix[j][k] += Jxw * basis_value[j][i] * basis_value[k][i];
            }
        }
	}
}

void Parabolic::buildMatrix_RHS_Sparse(Eigen::SparseMatrix<double> &A_matrix_sparse,
                                       Eigen::SparseMatrix<double> &M_matrix_sparse, 
                                       std::vector<double>& Rhs, double time)
{
    auto n_ele = mesh.n_element();
    auto n_pnt = mesh.n_point();
    A_matrix_sparse.resize(n_pnt,n_pnt);
    M_matrix_sparse.resize(n_pnt,n_pnt);
    Rhs.resize(n_pnt);
    // 定义一个用于存储稀疏矩阵的三元组
    std::vector<Eigen::Triplet<double>> Atriple, Mtriple;
    for (size_t i = 0; i != n_ele; ++i){
    	std::vector<int> NV = mesh.getEleVtx(i);
        std::vector<std::vector<double>> Ele_Amatrix(3);
        std::vector<std::vector<double>> Ele_Mmatrix(3);
        for (size_t j = 0; j != 3; ++j){
            Ele_Amatrix[j].resize(3,0.0);
            Ele_Mmatrix[j].resize(3,0.0);
        }
        std::vector<double> Ele_Rhs(3,0.0);
        buildEleMatrix_RHS(i,Ele_Amatrix,Ele_Mmatrix,Ele_Rhs,time);
        // 稀疏矩阵存储技术
        // 利用经典的三元组插入方式来存储A_matrix
        Atriple.push_back(Eigen::Triplet<double>(NV[0],NV[0],Ele_Amatrix[0][0]));
        Atriple.push_back(Eigen::Triplet<double>(NV[0],NV[1],Ele_Amatrix[0][1]));
        Atriple.push_back(Eigen::Triplet<double>(NV[0],NV[2],Ele_Amatrix[0][2]));
        Atriple.push_back(Eigen::Triplet<double>(NV[1],NV[0],Ele_Amatrix[1][0]));
        Atriple.push_back(Eigen::Triplet<double>(NV[1],NV[1],Ele_Amatrix[1][1]));
        Atriple.push_back(Eigen::Triplet<double>(NV[1],NV[2],Ele_Amatrix[1][2]));
        Atriple.push_back(Eigen::Triplet<double>(NV[2],NV[0],Ele_Amatrix[2][0]));
        Atriple.push_back(Eigen::Triplet<double>(NV[2],NV[1],Ele_Amatrix[2][1]));
        Atriple.push_back(Eigen::Triplet<double>(NV[2],NV[2],Ele_Amatrix[2][2]));
        // 利用经典的三元组插入方式来存储M_matrix
        Mtriple.push_back(Eigen::Triplet<double>(NV[0],NV[0],Ele_Mmatrix[0][0]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[0],NV[1],Ele_Mmatrix[0][1]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[0],NV[2],Ele_Mmatrix[0][2]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[1],NV[0],Ele_Mmatrix[1][0]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[1],NV[1],Ele_Mmatrix[1][1]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[1],NV[2],Ele_Mmatrix[1][2]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[2],NV[0],Ele_Mmatrix[2][0]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[2],NV[1],Ele_Mmatrix[2][1]));
        Mtriple.push_back(Eigen::Triplet<double>(NV[2],NV[2],Ele_Mmatrix[2][2]));

        // 不利用稀疏矩阵存储右端项Rhs
        Rhs[NV[0]] += Ele_Rhs[0];
        Rhs[NV[1]] += Ele_Rhs[1];
        Rhs[NV[2]] += Ele_Rhs[2];
    }
    // 把三元组转换成稀疏矩阵
    A_matrix_sparse.setFromTriplets(Atriple.begin(),Atriple.end());
    M_matrix_sparse.setFromTriplets(Mtriple.begin(),Mtriple.end());
}

// 对于稀疏矩阵A_matrix_sparse的边界条件处理(dirichlet边界条件)
void Parabolic::DealBoundary_Sparse(Eigen::SparseMatrix<double> &M_matrix_sparse, 
                                    std::vector<double>& Rhs, double t)
{
    auto n_pnt = mesh.n_point();
    auto n_bp = mesh.n_boundaryPoint();
    std::vector<int> BV = mesh.getBndPnt(); // 得到边界点
    for (size_t i = 0; i != n_bp; ++i){
        Point Pnt_i = mesh.getPnt(BV[i]); // 得到第i个单元的点
        double val = pde.u_boundary(Pnt_i,t);
        Rhs[BV[i]] = M_matrix_sparse.coeffRef(BV[i],BV[i]) * val;
        for (size_t j = 0; j != n_pnt; ++j){
            if (j != BV[i]){
                M_matrix_sparse.coeffRef(BV[i],j) = 0.0;
            }
        }
        for (size_t j = 0; j != n_pnt; ++j){
            if (j != BV[i]){
                Rhs[j] -= M_matrix_sparse.coeffRef(j, BV[i]) * val;
                M_matrix_sparse.coeffRef(j,BV[i]) = 0.0;
            }
        }
    }
    // 去除稀疏矩阵中的零元素
    M_matrix_sparse = M_matrix_sparse.pruned();
    // 压缩剩余空间
    M_matrix_sparse.makeCompressed();
}

// 全离散显格式
void Parabolic::fullDiscreteFormat(Eigen::SparseMatrix<double> &A_matrix_sparse,
                        const Eigen::SparseMatrix<double> &M_matrix_sparse,
                        const std::vector<double> &u_h, std::vector<double> &Rhs)
{
    // 去除矩阵A_matrix_sparse中的零元素
    A_matrix_sparse = A_matrix_sparse.pruned();
    // 压缩剩余空间
    A_matrix_sparse.makeCompressed();
    // 利用Eigen中的Matrix存储u_h,Rhs
    Eigen::MatrixXd u_h_matrix(u_h.size(),1);
    Eigen::MatrixXd Rhs_matrix(Rhs.size(),1);
    for (decltype(u_h.size()) i = 0; i != u_h.size(); ++i)
        u_h_matrix(i,0) = u_h[i];
    for (decltype(Rhs.size()) i = 0; i != Rhs.size(); ++i)
        Rhs_matrix(i,0) = Rhs[i];
    // Eigen中二元操作符支持稀疏矩阵和密集矩阵的混合操作
    // 实现全离散显格式
    Rhs_matrix = (M_matrix_sparse - Dt * A_matrix_sparse) * u_h_matrix + Rhs_matrix * Dt;
    // 把得到的Eigen密集型矩阵转换为vector类型
    for (decltype(Rhs.size()) i = 0; i != Rhs.size(); ++i)
        Rhs[i] = Rhs_matrix(i,0);
}

// 全离散隐格式
void Parabolic::fullDiscreteImplicitFormat(Eigen::SparseMatrix<double> &A_matrix_sparse,
                        Eigen::SparseMatrix<double> &M_matrix_sparse,
                        const std::vector<double> &u_h, std::vector<double> &Rhs)
{
    // 去除矩阵A_matrix_sparse中的零元素
    A_matrix_sparse = A_matrix_sparse.pruned();
    // 压缩剩余空间
    A_matrix_sparse.makeCompressed();
    // 利用Eigen中的Matrix存储u_h,Rhs
    Eigen::MatrixXd u_h_matrix(u_h.size(),1);
    Eigen::MatrixXd Rhs_matrix(Rhs.size(),1);
    for (decltype(u_h.size()) i = 0; i != u_h.size(); ++i)
        u_h_matrix(i,0) = u_h[i];
    for (decltype(Rhs.size()) i = 0; i != Rhs.size(); ++i)
        Rhs_matrix(i,0) = Rhs[i];
    // 实现全离散隐格式
    Rhs_matrix = Dt * Rhs_matrix + M_matrix_sparse * u_h_matrix;
    M_matrix_sparse = M_matrix_sparse + Dt * A_matrix_sparse;
    // 把得到的Eigen密集型矩阵转换为vector类型
    for (decltype(Rhs.size()) i = 0; i != Rhs.size(); ++i)
        Rhs[i] = Rhs_matrix(i,0);
}

double Parabolic::ComputerL2Error(const std::vector<double> &f, double time)
{
	double err=0.0;
	int n_ele=mesh.n_element();
	for(int j=0; j<n_ele; j++){
		std::vector<int> NV = mesh.getEleVtx(j);
		std::vector<Point> EP(3);
		EP[0]=mesh.getPnt(NV[0]);
		EP[1]=mesh.getPnt(NV[1]);
		EP[2]=mesh.getPnt(NV[2]);
		double volume = tmpEle.getArea();
		std::vector<Point> GP = tmpEle.getGaussPnt();

		std::vector<double> GW = tmpEle.getGaussWeight();
		std::vector<Point> q_point = tmpEle.Local_to_Global(GP, EP); 
		std::vector<double> jacobian = tmpEle.Local_to_Global_jacobian(GP, EP);
		std::vector<std::vector<double> > basis_value;
		basisValue(q_point, EP, basis_value);

		for(int i = 0; i < q_point.size(); i++){
			double Jxw = GW[i]*jacobian[i]*volume;
			double f_value = pde.u_exact(q_point[i],time);
			double u_h_val = f[NV[0]]*basis_value[0][i]+f[NV[1]]*basis_value[1][i]+f[NV[2]]*basis_value[2][i];
			double df_val = f_value-u_h_val;
			err += Jxw*df_val*df_val;
		}
	}
	err=sqrt(err);
	return err;
}

double Parabolic::ComputerH1Error(const std::vector<double> &f, double time)
{
   double err=0.0;
   int n_ele=mesh.n_element();
   for(int j=0; j<n_ele; j++){
        std::vector<int> NV=mesh.getEleVtx(j);
        std::vector<Point> EP(3);
        EP[0]=mesh.getPnt(NV[0]);
        EP[1]=mesh.getPnt(NV[1]);
        EP[2]=mesh.getPnt(NV[2]);
        std::vector<Point> GP=tmpEle.getGaussPnt();

        std::vector<double> GW=tmpEle.getGaussWeight();
        std::vector<Point> q_point = tmpEle.Local_to_Global(GP, EP); 
        std::vector<double> jacobian = tmpEle.Local_to_Global_jacobian(GP, EP);
        std::vector<std::vector<std::vector<double>>> basis_grad;
        basisGrad(q_point, EP, basis_grad);
	  
        for(int i = 0; i < q_point.size(); i++){
            double Jxw = GW[i]*jacobian[i];
            std::vector<double> f_value = pde.u_exact_grad(q_point[i], time);
            double u_h_grad1=0.0;
            double u_h_grad2=0.0;

            u_h_grad1 = f[NV[0]] * basis_grad[0][i][0] + f[NV[1]] * basis_grad[1][i][0] +
                f[NV[2]] * basis_grad[2][i][0];
            u_h_grad2 = f[NV[0]] * basis_grad[0][i][1] + f[NV[1]] * basis_grad[1][i][1] +
                f[NV[2]] * basis_grad[2][i][1];
            
            double df_val1 = f_value[0] - u_h_grad1;
            double df_val2 = f_value[1] - u_h_grad2;
            err += Jxw * df_val1 * df_val1;
            err += Jxw * df_val2 * df_val2;
        }
   }
   err = sqrt(err);
   return err;
}

// 输出u_h在边界处的值
void Parabolic::checkBoundary(const std::vector<double> &u_h, double t)
{
    std::vector<int> boundaryPointIndex = mesh.getBndPnt();
    for (size_t i = 0; i != mesh.n_boundaryPoint(); ++i){
        std::cout << u_h[boundaryPointIndex[i]] << "\t" << pde.u_boundary(mesh.getPnt(boundaryPointIndex[i]),t) << "\t"
                    << pde.u_boundary(mesh.getPnt(boundaryPointIndex[i]),t) - u_h[boundaryPointIndex[i]] << "\n";
    }
    std::cout << "\n";
}

void Parabolic::checkValue(const std::vector<double> &u_h, double t)
{
    for (size_t i = 0; i != mesh.n_point(); ++i){
        std::cout << u_h[i] << "\t" << pde.u_exact(mesh.getPnt(i), t)
                    <<"\t" << pde.u_exact(mesh.getPnt(i),t) - u_h[i] << "\n";
    }
    std::cout << std::endl;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////
// 类的接口函数
////////////////////////////////////////////////////////////////////////////////////////////////////////

void Parabolic::runFullDiscreteFormat(const std::string &outputFile)
{
    int n_pnt = mesh.n_point();             // 得到网格点的个数
    u_h.resize(n_pnt);
    // 给u_h赋初始值
    for (size_t i = 0; i != n_pnt; ++i)
        u_h[i] = pde.u_exact(mesh.getPnt(i),0);  
    double L2Error, H1Error;
    std::ofstream out(outputFile);
    out << "time\tL2Error\tH1Error" << "\n";  
    for (double time = Dt; time <= 1.0000009; time = time + Dt){
        // 定义一个用于存储A_matrix的稀疏矩阵A_matrix_sparse
        Eigen::SparseMatrix<double> A_matrix_sparse, M_matrix_sparse;
        std::vector<double> Rhs;
        // 求解刚度矩阵A(a(u,v)),M((u,v)),Rhs((f,v))
        buildMatrix_RHS_Sparse(A_matrix_sparse, M_matrix_sparse, Rhs, time);
        // 求解抛物方程的全离散显格式
        fullDiscreteFormat(A_matrix_sparse, M_matrix_sparse, u_h, Rhs);
        // 边界条件的处理
        DealBoundary_Sparse(M_matrix_sparse, Rhs, time);
        // k可以取1-5; gmres:case 3(求解PDE离散后形成的代数系统)
        std::cout << time << std::endl;
        SolveLinearASparse(M_matrix_sparse, u_h, Rhs, 3);
        L2Error = ComputerL2Error(u_h, time);      // 求解L2误差
        H1Error = ComputerH1Error(u_h, time);      // 求解H1误差 
        out << time << "\t" << L2Error << "\t" << H1Error << "\n";
    }
    out.close();
}

void Parabolic::runFullDiscreteImplicitFormat(const std::string &outputFile)
{
    int n_pnt = mesh.n_point();             // 得到网格点的个数
    u_h.resize(n_pnt);
    double L2Error, H1Error;
    // 给u_h赋初始值
    for (size_t i = 0; i != n_pnt; ++i)
        u_h[i] = pde.u_exact(mesh.getPnt(i),0); 
    L2Error = ComputerL2Error(u_h,0);
    H1Error = ComputerH1Error(u_h,0);
    std::ofstream out(outputFile);
    out << "time\tL2Error\tH1Error" << "\n";  
    for (double time = Dt; time <= 1.00009; time += Dt){
        // 定义一个用于存储A_matrix的稀疏矩阵A_matrix_sparse
        Eigen::SparseMatrix<double> A_matrix_sparse, M_matrix_sparse;
        std::vector<double> Rhs;
        // 求解刚度矩阵A(a(u,v)),M((u,v)),Rhs((f,v))
        buildMatrix_RHS_Sparse(A_matrix_sparse, M_matrix_sparse, Rhs, time);
        // 求解抛物方程的全离散隐格式
        fullDiscreteImplicitFormat(A_matrix_sparse, M_matrix_sparse, u_h, Rhs);
        // 边界条件的处理
        DealBoundary_Sparse(M_matrix_sparse, Rhs, time);
        // k可以取1-5; gmres:case 3(求解PDE离散后形成的代数系统)
        std::cout << time << std::endl;
        SolveLinearASparse(M_matrix_sparse, u_h, Rhs, 3);
        L2Error = ComputerL2Error(u_h, time);      // 求解L2误差
        H1Error = ComputerH1Error(u_h, time);      // 求解H1误差 
        out << time << "\t" << L2Error << "\t" << H1Error << "\n";
        //checkBoundary(u_h,time);
        //checkValue(u_h,time);
    }
    out.close();
}
<file_sep>/Possion方程-C++-p1元--Eigen(no sparse)/Matrix.h
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <fstream>
#include <Eigen/Dense>  // 线性代数库
#include <Eigen/IterativeLinearSolvers> // 迭代法求解线性方程组的求解库
#include <unsupported/Eigen/IterativeSolvers>

class Point
{
	friend std::ostream& operator<< (std::ostream&, const Point&);
	friend std::istream& operator>> (std::istream&, Point&);
	friend Point operator+ (const Point&, const Point&);
	friend Point operator- (const Point&, const Point&);
	friend Point midpoint(const Point&, const Point&);
	friend double distance(const Point&, const Point&);
    // 求三角形的质心,可以求加权质心
	friend Point barycenter(const std::vector<Point>&, const double *);
    friend double AREA(const Point&, const Point&, const Point&);

	private:
		double x[2];   
	
	public:
		Point();  
		Point(const double*);  
 		Point(const Point&);
 		Point(const double&, const double&);  

 		Point& operator= (const Point&);
 		double& operator[] (int);
 		const double& operator[] (int) const;
 		double length() const;
 		Point& operator+= (const Point&);		
 		Point& operator-= (const Point&);
 		Point& operator*= (const double&);
 		Point& operator/= (const double&);
};

class TmpEle
{
	friend std::ostream& operator<< (std::ostream&, const TmpEle&);

	private:
		std::vector<Point> Pnt;
		std::vector<Point> GaussPnt;
		std::vector<double> GaussWeight;

	public:
		// 构造函数
		TmpEle();
		TmpEle(int);
		void buildTE(int);

		Point getPnt(int);
		double getArea();
		std::vector<Point> getGaussPnt();
		std::vector<double> getGaussWeight();

		// 标准单元到普通单元
		Point Local_to_Global(const Point&, const std::vector<Point> &) const;
		// 普通单元到标准单元
		Point Global_to_Local(const Point&, const std::vector<Point> &) const;
		// 标准单元到普通单元
		std::vector<Point> Local_to_Global(const std::vector<Point> &,const std::vector<Point> & ) const;
		// 普通单元到标准单元
		std::vector<Point> Global_to_Local(const std::vector<Point>&, const std::vector<Point>&) const;
		// 从标准单元到普通单元的雅克比矩阵
		double Local_to_Global_jacobian(const Point&, const std::vector<Point>&) const;
		// 从普通单元到标准单元的雅克比矩阵
		double Global_to_Local_jacobian(const Point&, const std::vector<Point>&) const;
		// 从标准单元到普通单元的雅克比矩阵
		std::vector<double> Local_to_Global_jacobian(const std::vector<Point>&, const std::vector<Point>&) const;
		// 从普通单元到标准单元的雅克比矩阵
		std::vector<double> Global_to_Local_jacobian(const std::vector<Point>&, const std::vector<Point>&) const; 
};

class Mesh
{
	friend std::ostream& operator<< (std::ostream&, const Mesh&);

	private:
		std::vector<Point> Pnt;   // 三角单元的点
		std::vector<int> BndPnt;  // 三角单元的边界点
		std::vector<std::vector<int>> Ele;  // 三角单元的边

	public:
		Mesh() = default;  // 默认构造函数

		void readData(const std::string&);
		int getEleVtx(int , int);
		std::vector<int> getEleVtx(int);
		Point getPnt(int);
		std::vector<Point> getPnt(std::vector<int>&);
		std::vector<int> getBndPnt();
		int n_element(); // 得到三角单元的数量
		int n_point();  // 得到点的个数
		int n_boundaryPoint(); // 得到边界点的个数
};

class PDE
{	
	public:
		// 偏微分方程的边界条件
		double u_boundary(const Point&);

		// 偏微分方程的右端项
		double f(const Point&);

		// 偏微分方程解析解
		double u_exact(const Point&);
};


class Matrix
{
	public:
		// 数据结构
		PDE pde;	
		Mesh mesh;     
		TmpEle tmpEle;  
		std::vector<double> u_h;

		Matrix() = default;
		Matrix(const std::string&, int);  
		// gauss-seidel迭代
		void GaussSeidel(std::vector<std::vector<double>>&, std::vector<double>&, std::vector<double>&);
       
        // 利用Eigen库求解代数方程组
        void SolveLinearA(std::vector<std::vector<double>>&, std::vector<double>&, std::vector<double>&, int);

		// 基函数的构造
		void basisValue(Point&, std::vector<Point>&, std::vector<double>&);
		void basisValue(std::vector<Point>&, std::vector<Point>&, std::vector<std::vector<double>>&);

		// 梯度算子的计算
		void basisGrad(Point&, std::vector<Point>&, std::vector<std::vector<double>>&);
		void basisGrad(std::vector<Point>&, std::vector<Point>&, std::vector<std::vector<std::vector<double>>>&);

		// 单元刚度矩阵
		void buildEleMatrix_RHS(int, std::vector<std::vector<double>>&, std::vector<double>&);

		// 总刚度矩阵
		void buildMatrix_RHS(std::vector<std::vector<double>>&, std::vector<double>&);

		
		// 边界条件的处理
		void DealBoundary(std::vector<std::vector<double>>&, std::vector<double>&, std::vector<double>&);

		// 计算L2误差
		double ComputerL2Error(std::vector<double> &);
	
};
<file_sep>/Heat方程-C++显格式/PDE.cpp
/*************************************************************************
	> File Name: PDE.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时45分37秒
 ************************************************************************/

#include "PDE.h"

// 边界条件(第一类边界条件)
double PDE::u_boundary(const Point &p, double t)
{
	double value;
    value = 0.0;
    return value;
}

double PDE::f(const Point &p, double t)
{
	double value;
    value = exp(t) * ( (p[0] * p[0] - p[0]) * (p[1] * p[1] - p[1]) - 
                    2 * (p[0] * p[0] - p[0] + p[1] * p[1] - p[1]) ); 
	return value; 
}

double PDE::u_exact(const Point& p, double t)
{
	double value;
	value = (p[0] * p[0] - p[0]) * (p[1] * p[1] - p[1]) * exp(t);
    return value;
}

std::vector<double> PDE::u_exact_grad(const Point &p, double t)
{
    std::vector<double> val(2);
    val[0] = -exp(t) * (-p[1] * p[1] + p[1]) * (2 * p[0] - 1);
    val[1] = -exp(t) * (-p[0] * p[0] + p[0]) * (2 * p[1] - 1);
    return val;
}
<file_sep>/Possion方程--C++(CVT自适应)/Makefile
all: CVT_Adaptive_finite

CVT_Adaptive_finite: CVT_Adaptive_finite.o Possion.o Matrix.o density.o
	g++ -std=c++0x -o CVT_Adaptive_finite CVT_Adaptive_finite.o Possion.o Matrix.o density.o

CVT_Adaptive_finite.o: CVT_Adaptive_finite.cpp Possion.h density.h
	g++ -std=c++0x -c CVT_Adaptive_finite.cpp

Possion.o: Possion.cpp Possion.h
	g++ -std=c++0x -c Possion.cpp

Matrix.o: Matrix.cpp Matrix.h
	g++ -std=c++0x -c Matrix.cpp

density.o: density.cpp density.h
	g++ -std=c++0x -c density.cpp
<file_sep>/Heat方程-C++显格式/Point.h
/*************************************************************************
	> File Name: Point.h
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 08时51分59秒
 ************************************************************************/

#ifndef POINT_H
#define POINT_H

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>  // 线性代数库
#include <Eigen/Sparse> // 稀疏存储矩阵库,以及稀疏存储的线性方程组求解库
#include <Eigen/IterativeLinearSolvers> // 迭代法求解线性方程组的求解库
#include <unsupported/Eigen/IterativeSolvers>

class Point
{
	friend std::ostream& operator<< (std::ostream&, const Point&);
	friend std::istream& operator>> (std::istream&, Point&);
	friend Point operator+ (const Point&, const Point&);
	friend Point operator- (const Point&, const Point&);
	friend Point midpoint(const Point&, const Point&);
	friend double distance(const Point&, const Point&);
    // 求三角形的质心,可以求加权质心
	friend Point barycenter(const std::vector<Point>&, const double *);
    friend double AREA(const Point&, const Point&, const Point&);

private:
	double x[2];   
	
public:
	Point();  
	Point(const double*);  
 	Point(const Point&);
 	Point(const double&, const double&);  
 	Point& operator= (const Point&);
 	double& operator[] (int);
 	const double& operator[] (int) const;
 	double length() const;
 	Point& operator+= (const Point&);		
 	Point& operator-= (const Point&);
 	Point& operator*= (const double&);
 	Point& operator/= (const double&);
};

#endif
<file_sep>/Possion方程-C++-p1元--Eigen(no sparse)/main.cpp
#include "Matrix.h"

using std::cout;
using std::endl;
using std::cin;

int main()
{
	Matrix Possion("Trimesh.txt",7);
	PDE pde;

	int n_pnt = Possion.mesh.n_point();  // 得到网格点的个数
	Possion.u_h.resize(n_pnt);         

	std::vector<std::vector<double>> A_matrix;
	std::vector<double> Rhs;

	Possion.buildMatrix_RHS(A_matrix,Rhs);   // 得到总的左边刚度矩阵和右端刚度矩阵

	Possion.DealBoundary(A_matrix, Possion.u_h, Rhs);  // 边界处理

	//Possion.GaussSeidel(A_matrix,Possion.u_h, Rhs);  // 高斯赛德尔迭代
    // k可以取1到14,取整数,每一个数字对应一种计算方法,
    // 1--9 为直接解法; 10--14为迭代解法;GMRES,case12 优先选择
    Possion.SolveLinearA(A_matrix,Possion.u_h,Rhs,12);


	double L2Error=Possion.ComputerL2Error(Possion.u_h);  // 求解L2误差
   	std::cout<<"L2Error="<<L2Error<<std::endl;

   	// 把数值解u_h和误差L2写入到文件中
   	std::ofstream out("Results.txt");
   	for (auto i : Possion.u_h)
   		out << i << "\n";
   	out.close();

   	return 0;
}
<file_sep>/Possion方程-C++-P1元--Eigen(Sparse)/main.cpp
#include "Matrix.h"

using std::cout;
using std::endl;
using std::cin;

int main()
{
	Matrix Possion("Trimesh.txt",7);
	PDE pde;

	int n_pnt = Possion.mesh.n_point();  // 得到网格点的个数
    Possion.u_h.resize(n_pnt);

	std::vector<std::vector<double>> A_matrix;
	std::vector<double> Rhs;

	//Possion.buildMatrix_RHS(A_matrix,Rhs);   // 得到总的左边刚度矩阵和右端刚度矩阵
   
    // 定义一个用于存储A_matrix的稀疏矩阵A_matrix_sparse
    Eigen::SparseMatrix<double> A_matrix_sparse(A_matrix.size(),A_matrix.size());
    Possion.buildMatrix_RHS_Sparse(A_matrix_sparse,Rhs);
    
    Possion.DealBoundary_Sparse(A_matrix_sparse, Possion.u_h, Rhs);
        
    // k 可以取1-5 ; gmres:case 3
    Possion.SolveLinearASparse(A_matrix_sparse,Possion.u_h,Rhs,3);

    double L2Error=Possion.ComputerL2Error(Possion.u_h);  // 求解L2误差
   	std::cout<<"L2Error="<<L2Error<<std::endl;

    // 把数值解u_h和误差L2写入到文件中
   	std::ofstream out("Results.txt");
   	for (auto i : Possion.u_h)
   		out << i << "\n";
   	out.close();

    return 0;
}
<file_sep>/Heat方程-C++隐格式-非零边界/PDE.cpp
/*************************************************************************
	> File Name: PDE.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 09时45分37秒
 ************************************************************************/

#include "PDE.h"

// 边界条件(第一类边界条件)
double PDE::u_boundary(const Point &p, double t)
{
	double value;
    value = exp(0.5 * (p[0] + p[1]) - t);
    return value;
}

double PDE::f(const Point &p, double t)
{
	double value;
	value = -3/2.0 * exp(0.5 * (p[0] + p[1]) - t);
    return value; 
}

double PDE::u_exact(const Point& p, double t)
{
	double value;
    value = exp(0.5 * (p[0] + p[1]) - t);
    return value;
}

std::vector<double> PDE::u_exact_grad(const Point &p, double t)
{
    std::vector<double> val(2);
    val[0] = 0.5 * exp(0.5 * (p[0] + p[1]) - t);
    val[1] = 0.5 * exp(0.5 * (p[0] + p[1]) - t);
    return val;
}
<file_sep>/Heat方程-C++显格式/Test.cpp
/*************************************************************************
	> File Name: Test.cpp
	> Author: ZhirunZheng 
	> Mail: <EMAIL> 
	> Created Time: 2018年10月22日 星期一 15时07分53秒
 ************************************************************************/
#include "Parabolic.h"

int main()
{
    int index = 0;
    for (int i = 0; i != 5; ++i){
        std::cout << "*********************"<<i<<"******************************"<<std::endl;
        std::string inputFile = "Trimesh_";
        std::string outputFile = "Error_";
        inputFile = inputFile + std::to_string(index + i) + ".txt";
        outputFile = outputFile + std::to_string(index + i) + ".txt";
        Parabolic solve(inputFile, 5);
        solve.runFullDiscreteFormat(outputFile);        // 全离散显格式
        Dt = Dt / 4.0;      // 二维情况下的网格比要小于1/4
    }
    return 0;
}
