MoM矩量法(一):剖分数据的处理

MoM矩量法(一):剖分数据的处理

写在开头:沉寂了一年多,终于又要开始写博客了。这一年多来经历了很多事情,我也成长了不少,从Java到机器学习再到计算电磁学,方向的转变还是挺大的。虽然方向在变,但之前写代码留下的功底作用还是非常大的。好了,废话不多说,开始进入主题吧。
矩量法的计算理论还有公式推导很多专门的书上都有,这里不再赘述。需要注意的是,要看懂书上的推导过程还是需要比较深的物理和数学功底的,所以一开始看的比较吃力的话,不要担心,慢慢来。
这里采用的公式来自于聂在平主编的《目标与环境 电磁散射特性建模——理论、方法与实现(应用篇)》,所有数据均采用国际标准单位。公式的推导就不再多说了,直接开始coding吧!

完整代码已贴在文章最后。文章和代码思路借鉴了https://blog.csdn.net/u014411646/article/details/99689460,谢谢原作者!

一、剖分数据

初始剖分数据包括两个文件,分别是包含点元素信息的文件pdata.txt和包含三角单元信息的文件tdata.txt,先一起来看一看初始数据。

1.1点元素数据文件pdata.txt

在这里插入图片描述
横向为坐标值,纵向为点元素序号。例如,上图中显示第4个点的坐标为(-3.7033783e-02, 6.4144393e-02, 4.4386259e-01)。所以pdata.txt的行个数就代表点元素的个数。

1.2三角单元数据文件tdata.txt

在这里插入图片描述
横向为点元素的序号,纵向为三角单元序号。例如,从上图中可以知道第3个三角单元的三个顶点分别为1,4和5。则tdata.txt的行个数即为三角单元的个数。

上面两组数据的意义是非常明显的,例如,从第二组数据中我们可以得到第2个三角单元的三个顶点分别是1,3和4,而从第一组数据又可以得到第1号点的坐标为(0.0000000e+00 0.0000000e+00 4.5000000e-01),第3号点的坐标为(3.7033783e-02 6.4144393e-02 4.4386259e-01),第4号点的坐标为(-3.7033783e-02 6.4144393e-02 4.4386259e-01)。所以,第2个三角单元三个顶点坐标分别为(0.0000000e+00 0.0000000e+00 4.5000000e-01)、(3.7033783e-02 6.4144393e-02 4.4386259e-01)和(-3.7033783e-02 6.4144393e-02 4.4386259e-01)。

二、初始数据的处理

初始数据的处理其实有很多种方法,但核心都是为了之后计算阻抗矩和激励矩阵的简便。
接下来对一些核心代码进行讲解:

mat Edge_list;              //[端点1, 端点2, 对点1, 对点2, 左三角单元, 右三角单元]
mat TriCenter;
vec EdgeLength;
vec TriArea;
int TriTotal, EdgeTotal;
mat pdata, tdata;
  1. TriTotal为三角单元的个数,EdgeTotal为边的个数;
  2. Edge_list是一个维度为(EdgeTotal,6)的二维矩阵,其中储存的是剖分数据中边的信息,行代表第几条边,列从左到右依次代表[第一个端点, 第二个端点, 与边相对的第一个端点,与边相对的第二个端点, 左边的三角单元, 右边的三角单元];
  3. TriCenter是(TriTotal,3)的二维矩阵,储存的是三角单元的重心坐标;
  4. EdgeLength是长度为EdgeTotal向量,储存的是对应边的长度;
  5. TriArea是长度为TriTotal向量,储存的是对应三角单元的面积;
  6. pdata, tdata为原始的pdata.txt和tdata.txt初始数据。
下面贴上代码:
/***********************************************************/
/*         该函数主要用来对输入的初始数据进行处理            */
/*                                                         */
/*                                                         */
/***********************************************************/
/*                                                         */
/*                 Author:Tianshaowenwen                  */
/*                                                         */
/***********************************************************/
#include "rwg1.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "armadillo"

using namespace arma;
using namespace std;


extern string p_filename, t_filename;


mat Edge_list;              //[端点1, 端点2, 对点1, 对点2, 左三角单元, 右三角单元]
mat TriCenter;
vec EdgeLength;
vec TriArea;
int TriTotal, EdgeTotal;
mat pdata, tdata;

void rwg1_main(){
	pdata.load(p_filename, raw_ascii);
	tdata.load(t_filename, raw_ascii);
	TriTotal = tdata.n_rows;
	edgeNumber();
	edgeElement();
	area();
	Edge_list.save("Edge_list.txt", raw_ascii);
	TriCenter.save("TriCenter.txt", raw_ascii);
	EdgeLength.save("EdgeLength.txt", raw_ascii);
	TriArea.save("TriArea.txt", raw_ascii);
	cout << "\n" << "TriTotal= " << TriTotal << endl;
	cout << "\n" << "EdgeTotal= " << EdgeTotal << endl;
}
/***********************************************************/
/*              Function:void edgeNumber()                */
/***********************************************************/
/* 1.统计边的个数                                          */
/***********************************************************/
void edgeNumber(){
	int edgetotal = 0;
	for (int i = 0;i < TriTotal;i++) {
		int N[3] = { 0 }, M[3] = { 0 };
		N[0] = tdata(i, 0) - 1;
		N[1] = tdata(i, 1) - 1;
		N[2] = tdata(i, 2) - 1;
		for (int j = i + 1;j < TriTotal;j++) {
			M[0] = tdata(j, 0) - 1;
			M[1] = tdata(j, 1) - 1;
			M[2] = tdata(j, 2) - 1;
			int sameNode = 0;
			for(int i1 = 0;i1 < 3;i1++) {
				for (int i2 = 0;i2 < 3;i2++) {
					if (N[i1] == M[i2]) {
						sameNode += 1;
					}
				}
			}
			if (sameNode == 2) {
				edgetotal += 1;
			}
		}
	}
	EdgeTotal = edgetotal;
}

/***********************************************************/
/*             Function:void edgeelement()                */
/***********************************************************/
/* 1.绑定edge两端节点和相对节点和左右三角单元                */
/* 2.求出edge长度                                           */
/*                                                         */
/***********************************************************/
void edgeElement(){
	vec edgeLength(EdgeTotal);
	mat edge_list(EdgeTotal, 6);
	int edgetotal = 0;
	for (int i = 0;i < TriTotal;i++) {
		int N[3] = { 0 }, M[3] = { 0 };
		N[0] = tdata(i, 0) - 1;
		N[1] = tdata(i, 1) - 1;
		N[2] = tdata(i, 2) - 1;
		for (int j = i + 1;j < TriTotal;j++) {
			M[0] = tdata(j, 0) - 1;
			M[1] = tdata(j, 1) - 1;
			M[2] = tdata(j, 2) - 1;
			int sameNode = 0;
			int same[2] = { 0 };
			for (int i1 = 0;i1 < 3;i1++) {
				for (int i2 = 0;i2 < 3;i2++) {
					if (N[i1] == M[i2]) {
						same[sameNode] = N[i1];
						sameNode += 1;
					}
				}
			}
			if (sameNode == 2) {
				edge_list(edgetotal, 0) = same[0];
				edge_list(edgetotal, 1) = same[1];
				for (int j = 0;j < 3;j++) {
					if (N[j] != same[0] && N[j] != same[1]) {
						edge_list(edgetotal, 2) = N[j];
					}
					if (M[j] != same[0] && M[j] != same[1]) {
						edge_list(edgetotal, 3) = M[j];
					}
				}
				edge_list(edgetotal, 4) = i;
				edge_list(edgetotal, 5) = j;

				double len = getDist(pdata.row(same[0]), pdata.row(same[1]));
				edgeLength(edgetotal) = len;
				edgetotal += 1;
			}
		}
	}
	Edge_list = edge_list;
	EdgeLength = edgeLength;
}
/***********************************************************/
/*                 Function:void area()                  */
/***********************************************************/
/* 1.计算三角单元的面积                                     */
/* 2.计算三角单元的重心                                     */
/***********************************************************/
void area() {
	vec triArea(TriTotal);
	mat triCenter(TriTotal, 3);
	float test = 0, test1 = 0;
	vec vector1(3), vector2(3);
	int N[3] = { 0 };
	for (int ia = 0;ia < TriTotal;ia++)
	{
		for (int i = 0;i < 3;i++)
			N[i] = tdata(ia, i) - 1;
		for (int j = 0;j < 3;j++)
		{
			vector1(j) = pdata(N[0], j) - pdata(N[1], j);
			vector2(j) = pdata(N[2], j) - pdata(N[1], j);
		}
		vec vector3 = cross(vector1, vector2);
		triArea(ia) = 0.5 * sqrt(vector3(0) * vector3(0) + vector3(1) * vector3(1) + vector3(2) * vector3(2));
		triCenter(ia, 0) = (1.0 / 3.0) * (pdata(N[0], 0) + pdata(N[1], 0) + pdata(N[2], 0));
		triCenter(ia, 1) = (1.0 / 3.0) * (pdata(N[0], 1) + pdata(N[1], 1) + pdata(N[2], 1));
		triCenter(ia, 2) = (1.0 / 3.0) * (pdata(N[0], 2) + pdata(N[1], 2) + pdata(N[2], 2));
	}
	TriArea = triArea;
	TriCenter = triCenter;
}
/***********************************************************/
/*                 Function:double getDist()              */
/***********************************************************/
/* 1.求两个向量的距离                                       */
/***********************************************************/
double getDist(rowvec r1, rowvec r2) {
	double dx = r1(0) - r2(0);
	double dy = r1(1) - r2(1);
	double dz = r1(2) - r2(2);
	double len = sqrt(dx * dx + dy * dy + dz * dz);
	return len;
}
/***********************************************************/
/*                 Function:rowvec cross()                */
/***********************************************************/
/* 1.求两个向量的差乘                                       */
/***********************************************************/
rowvec cross(rowvec r1,rowvec r2){
	rowvec res(3);
	res(0) = r1(1) * r2(2) - r1(2) * r2(1);
	res(1) = r1(2) * r2(0) - r1(0) * r2(2);
	res(2) = r1(0) * r2(1) - r1(1) * r2(0);
	return res;
}
/***********************************************************/
/*                 Function:double getDot()               */
/***********************************************************/
/* 1.求两个向量的点乘                                       */
/***********************************************************/
double getDot(rowvec r1, rowvec r2) {
	return r1(0) * r2(0) + r1(1) * r2(1) + r1(2) * r2(2);
}
/***********************************************************/
/*                 Function:vec getUnit()                 */
/***********************************************************/
/* 1.求单位向量                                            */
/***********************************************************/
rowvec getUnit(rowvec r1) {
	rowvec res(3);
	res(0) = 0.0;
	res(1) = 0.0;
	res(2) = 0.0;
	double len = sqrt(r1(0) * r1(0) + r1(1) * r1(1) + r1(2) * r1(2));
	if (len != 0.0) {
		res(0) = r1(0) / len;
		res(1) = r1(1) / len;
		res(2) = r1(2) / len;
	}
	return res;
}

今天写的这第一部分其实就是对初始数据的一些处理和需要用到的工具函数的创建。矩量法的重头其实还是阻抗矩阵的填充和矩阵方程的求解,关于阻抗矩阵的填充我会放在下次讲。

  • 7
    点赞
  • 19
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
以下是用矩量在MATLAB中剖分一个圆锥的代码: ``` % 定义圆锥的参数 r = 5; % 圆锥底面半径 h = 10; % 圆锥高度 % 定义剖分的参数 N = 30; % 切割成N层 M = 60; % 每层切割成M份 % 生成圆锥顶点 vertices = zeros(N+2, 3); vertices(1,:) = [0, 0, h]; % 生成圆锥底面的顶点 theta = linspace(0, 2*pi, M+1)'; theta(end) = []; x = r * cos(theta); y = r * sin(theta); z = zeros(size(theta)); vertices(end,:) = [0, 0, 0]; vertices(2:end-1,:) = [x, y, z]; % 生成三角形面片 faces = zeros(N*M, 3); f = 1; for i = 1:N for j = 1:M if j == M faces(f,:) = [M*(i-1)+j, M*(i-1)+1, M*i+1]; else faces(f,:) = [M*(i-1)+j, M*(i-1)+j+1, M*i+j+1]; end f = f + 1; end end % 绘制圆锥 figure; patch('Vertices',vertices,'Faces',faces,'FaceColor','g','EdgeColor','black'); axis equal; ``` 解释: - 首先,我们定义了圆锥的参数:半径为$r$,高度为$h$。 - 然后,我们定义了剖分的参数:将圆锥切割成$N$层,每层再切割成$M$份。 - 接下来,我们生成圆锥顶点和底面顶点,并存储在$vertices$矩阵中。$vertices$矩阵的第一行为圆锥顶点,第二行到第$N+1$行为底面顶点,第$N+2$行为底面中心。 - 最后,我们根据顶点生成三角形面片,并存储在$faces$矩阵中。$faces$矩阵的每一行包含三个顶点的索引,表示一个三角形面片。 - 绘制圆锥:使用MATLAB的patch函数,将$vertices$矩阵和$faces$矩阵作为输入参数,绘制出圆锥。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值