本文介绍了使用C++编程解决二维任意区域（含洞）内的拉普拉斯方程数值解。程序设计包括读取文本文件中的几何区域和边界信息，生成网格，组装刚度矩阵，并给出利用Matlab绘制结果的方案。文章强调了软件设计和接口的重要性，以及对有限元方法的理解。 摘要生成于 ，由 DeepSeek-R1 满血版支持， \nabla^2u(x, y) = 0 \ \ \ \ in \ \ \Omega \ \ \ (1) \\ u = \bar{u}(x,y) \ \ \ \ on \ \ \ \Gamma_u \ \ \ (2) \\ q = \frac{\partial u}{\partial \vec{n}} = 0 \ \ \ on \ \ \ \Gamma_q \ \ \ (3) \\ ∇ 2 u ( x , y ) = 0 in Ω ( 1 ) u = u ˉ ( x , y ) o n Γ u ​ ( 2 ) q = ∂ n ∂ u ​ = 0 o n Γ q ​ ( 3 )

// 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 ; 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剖二维网格教程附代码 。