有限元方法求解二维拉普拉斯方程C++实现

前言

本文利用C++语言实现在二维任意区域(内部可有“洞)求解拉普拉斯方程的数值解。可以将区域以及边界信息写在文本文件中， 程序有专门的接口读入数据。

问题

区域

方程

$$\begin{gathered} {\nabla }^{2}u(x,y)=0in\Omega (1) \\ u=\bar{u}(x,y)on{\Gamma }_u(2) \\ q=\frac{\partial u}{\partial \vec{n}}=0on{\Gamma }_q(3) \end{gathered}$$

程序设计

几何区域

//
//  Created by guihun on 2022/10/24.
//

#ifndef Geometry_h
#define Geometry_h


#include <iostream>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <math.h>
class Point
{
/**
  * Point coordinate information
  */
    public:
        Point(const double x = 0, const double y = 0, const double z = 0):_x(x), _y(y) {}
        Point(const Point& p) : _x( p.getX() ), _y( p.getY() ) {;}
        ~Point() {}
        const double getX() const {return _x;}
        const double getY() const {return _y;}
        void setX(const double x) {_x = x;}
        void setY(const double y) {_y = y;}
        bool operator < (const Point& p) const;
        bool operator == (const Point& p) const;
        
    private:
        double _x, _y;
};


class Polygon
{
    public:
        Polygon() : _name(""), _vertex() { }
        Polygon(const string& name, const vector<Point>& vertex) : _name(name), _vertex(vertex) {}
        ~Polygon();
        const vector<Point>& getVertex() const {return _vertex;}
        const string& getName() const {return _name;}
        void clear() {_vertex.clear();}
        void append(const Point& p) {_vertex.push_back(p);}
        bool setPoint(const unsigned int index, const Point& p);
        const double getArea() const;
        bool isInPolygon(const Point& p) const;
		bool isOnPolygonEdge(const Point& p) const;
  
    private:
        string _name;
        vector<Point> _vertex;
};



class Region
{
    public:
        Region() : _polys(),
                   _window(){}
        ~Region();
    public:
        bool initByFile(const string& inputFile);

    public:
        const vector<Polygon>& getPolys() const {return _polys;}
        const Polygon& getWindow() const {return _window;}

    private:
    	// the hole of region
        vector<Polygon> _polys;
        Polygon _window;

};
#endif /* Geometry_h */

网格单元

//
//  Created by guihun on 2022/10/24.
//

#ifndef FEMMesh_h
#define FEMMesh_h

#include <iostream>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <math.h>

#include "Geometry.h"

using namespace std;

class Triangle
{

    /**
     * Triangle vertex index
     * _id[0] ~ _id[2] : three vertices index of triangle
     * _id[3] is the midpoint between _id[0] and _id[1]
     * _id[4] is the midpoint between _id[1] and _id[2]
     * _id[5] is the midpoint between _id[2] and _id[0]
     */
    public:
        Triangle() : _ids() {}
        explicit Triangle(const vector<unsigned int>& ids) : _ids(ids) {}
        ~Triangle() {}

        const vector<unsigned int>& getIds() const {return _ids;}
        void append(unsigned int id) {_ids.push_back(id);}
        void resize(const unsigned int size) {_ids.resize(size);}
        void clear() {_ids.clear();}
        bool setId(const unsigned int i, const unsigned int id);


    private:
        vector<unsigned int> _ids;
};


class IdTriangle : public Triangle
{
		/**
		 * TriangleId, TriangleNodesId, (the first two points(_id[0], _id[1] ->) are on the conductor boundary)
		 * _id[0], _id[1] storage 0 or 1 or 2
		 * _id[2] storage 3 or 4 or 5
		 * Example:
		 *      if _id[0] = 0 and _id[1] = 2 we can get _id[2] = 5
		 *      That is if we know the information stored in _id[0] and _id[1], 
		 *      we must know the information stored in _id[2]
		 *
		 */
		public:
				IdTriangle() :_triId(-1), _normal() {}
				IdTriangle(const vector<unsigned int>& ids, const unsigned int triId, const Point& normal): Triangle(ids), _triId(triId), _normal(normal) {}
				~IdTriangle() {}
				const unsigned int getTrId() const {return _triId;}
				void setTriId(const unsigned int triId) {_triId = triId;}
				const Point& getNormal() const {return _normal;}
				void setNormal(const Point& normal) {_normal = normal;}
    private:
      unsigned int _triId;
      Point _normal;

};


class FEMMesh
{
		/**
		 * mesh information of region 
		 */
		public:
			FEMMesh() : _nodes(), _elements(), _holerNodesIds(), _idTris(){}  
			~FEMMesh() {}
			const vector<Point>& getNodes() const {return _nodes;}
			const vector<Triangle>& getElements() const {return _elements;}			
			bool generateTriangleMesh(const Region& region);
			// convert triangle(3 nodes for p1 element) to triangle(6 nodes for p2 element)
			bool convertTriangle();
		public:
			// write mesh information to txt file;
			bool writeMesh(const string& location, const string& fileName) const;

		private:
			// all mesh nodes
			vector<Point> _nodes;
			// all triangles nodes index
			vector<Triangle> _elements;
			/**
		  	  * holes nodes index
			  * key : hole's idName
			  * value : node index which node on the hole edge.
			  */
			map<string, set<unsigned int> > _holerNodesIds;

			/**
			  * holes surround triangle
			  * key : hole's idName
			  * value : holes surround triangle information
			  */
			map<string, vector<IdTriangle> > _idTris;

};




#endif /* FEMMesh_h */

刚度矩阵组装

//
//  Created by guihun on 2022/10/24.
//
#ifndef _FEM_EQUATION_H_
#define _FEM_EQUATION_H_


#include <iostream>
#include <vector>
#include <map>

#include "FEMMesh.h"

using namespace std;

class FEMEquation
{
		public:
//.....

		private:
			void initGaussPara();

		private:
			unsigned int _dim;
			// Stiffness matrix
			map<unsigned int, map<unsigned int, double> > _stifMat;  
			// Load vector
			map<unsigned int, double> _loadMap;
			// solution
			vector<double> _soluVec;
			vector<double> _weight;
			vector<vector<double> > _gaussPoint;
};

#endif

数值结果

问题区域网格

u 值图像(结果导出借助Matlab画图)

总结

利用C++实现可以很好的学习软件设计思路，做各式适合接口(例如将几何区域的信息写入文本文件， 通过程序读取)，同时对有限元方法的单刚组总刚有一个更深的理解。本程序借助gmsh的C++接口实现剖网格， 可以参考我的另一篇博客 Gmsh剖二维网格教程附代码 。