有限元方法求解二维拉普拉斯方程C++实现_外汇行情

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

迪丽瓦拉

2025-05-31 16:05:22

0次

文章目录

- 前言
- 问题
- - 区域
  - 方程
- 程序设计
- - 几何区域
  - 网格单元
  - 刚度矩阵组装
  - 数值结果
  - - 问题区域网格
    - u 值图像(结果导出借助Matlab画图)
  - 总结

前言

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

问题

区域

在这里插入图片描述

方程

∇2u(x,y)=0inΩ(1)u=uˉ(x,y)onΓu(2)q=∂u∂n⃗=0onΓq(3)\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) \\ ∇2u(x,y)=0 in Ω (1)u=uˉ(x,y) on Γu (2)q=∂n∂u=0 on Γq (3)

程序设计

几何区域

//
//  Created by guihun on 2022/10/24.
//#ifndef Geometry_h
#define Geometry_h#include 
#include 
#include 
#include 
#include 
#include 
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& vertex) : _name(name), _vertex(vertex) {}~Polygon();const vector& 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 _vertex;
};class Region
{public:Region() : _polys(),_window(){}~Region();public:bool initByFile(const string& inputFile);public:const vector& getPolys() const {return _polys;}const Polygon& getWindow() const {return _window;}private:// the hole of regionvector _polys;Polygon _window;};
#endif /* Geometry_h */

网格单元

//
//  Created by guihun on 2022/10/24.
//#ifndef FEMMesh_h
#define FEMMesh_h#include 
#include 
#include 
#include 
#include 
#include #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& ids) : _ids(ids) {}~Triangle() {}const vector& 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 _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& 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& getNodes() const {return _nodes;}const vector& 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 nodesvector _nodes;// all triangles nodes indexvector _elements;/*** holes nodes index* key : hole's idName* value : node index which node on the hole edge.*/map > _holerNodesIds;/*** holes surround triangle* key : hole's idName* value : holes surround triangle information*/map > _idTris;};#endif /* FEMMesh_h */

刚度矩阵组装

//
//  Created by guihun on 2022/10/24.
//
#ifndef _FEM_EQUATION_H_
#define _FEM_EQUATION_H_#include 
#include 
#include #include "FEMMesh.h"using namespace std;class FEMEquation
{public:
//.....private:void initGaussPara();private:unsigned int _dim;// Stiffness matrixmap > _stifMat;  // Load vectormap _loadMap;// solutionvector _soluVec;vector _weight;vector > _gaussPoint;
};#endif