gctl_toolkits/archive/mshinterpolate/mshInterpolate.h

#ifndef _MSHINTERPOLATE_H
#define _MSHINTERPOLATE_H
#include "dataStruct.h"

class mshInterpolate
{
public:
	mshInterpolate(){}
	~mshInterpolate(){}
	int routine();
	int readMsh(char*,char*);
	int readList(char*);
	int initRef(char*);
	void interpolate_topo();
	void interpolate_linear();
	int outList();
private:
	vertexArray mshVert;
	triangleArray mshTri;
	_1dArray vertVal;

	vertexArray listVert;
	_1iArray listIds;
	_1dArray listVal;

	double refr,refR;
};

int mshInterpolate::readMsh(char* filename,char* valName)
{
	int temp_int,ele_type,attri_num,temp_attri,temp_id;
	double temp_val;
	vertex temp_vert;
	triangle temp_tri;
	string temp_str;
	stringstream temp_ss;

	ifstream mshin;
	if (open_infile(mshin,filename)) return -1;

	while(getline(mshin,temp_str))
	{
		if (temp_str == "$Nodes")
		{
			getline(mshin,temp_str);
			temp_ss.str("");
			temp_ss.clear();
			temp_ss << temp_str;
			temp_ss >> temp_int;
			for (int i = 0; i < temp_int; i++)
			{
				getline(mshin,temp_str);
				temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
				temp_ss >> temp_vert.id >> temp_vert.posic.x >> temp_vert.posic.y >> temp_vert.posic.z;
				temp_vert.set(temp_vert.posic);
				mshVert.push_back(temp_vert);
			}
			//顶点读入完成 初始化数组
			vertVal.resize(mshVert.size(),0.0);
		}
		else if (temp_str == "$Elements")
		{
			getline(mshin,temp_str);
			temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
			temp_ss >> temp_int;
			for (int i = 0; i < temp_int; i++)
			{
				getline(mshin,temp_str);
				temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
				temp_ss >> temp_tri.id >> ele_type; //这里temp_int是元素的id
				if (ele_type == 2) //只读入三角形
				{
					temp_ss >> attri_num;
					//第一个属性为物理组
					temp_ss >> temp_tri.phys;
					for (int a = 0; a < attri_num-1; a++) //忽略后面的属性值
						temp_ss >> temp_attri;
					//读入三棱锥顶面三角形顶点索引
					temp_ss >> temp_tri.vec[0]
							>> temp_tri.vec[1]
							>> temp_tri.vec[2];
					mshTri.push_back(temp_tri);
				}
			}
		}
		else if (temp_str == "$NodeData")
		{
			for (int i = 0; i < 2; i++) //先读入元素块的名称 按照关键字解析不同类型的元素值
				getline(mshin,temp_str);

			if (!strcmp(temp_str.c_str(),valName))
			{
				for (int i = 0; i < 6; i++) //跳过元素属性前面的值 最后一次为当前元素块的个数
				getline(mshin,temp_str);
				temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
				temp_ss >> temp_int;
				for (int i = 0; i < temp_int; i++)
				{
					getline(mshin,temp_str);
					temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
					temp_ss >> temp_id >> temp_val;
					vertVal[temp_id] = temp_val;
				}
			}
		}
		else continue;
	}
	mshin.close();
	return 0;
}

int mshInterpolate::readList(char* filename)
{
	int temp_id;
	vertex temp_vert;
	string temp_str;
	stringstream temp_ss;

	ifstream infile;
	if (open_infile(infile,filename)) return -1;

	while(getline(infile,temp_str))
	{
		if (*(temp_str.begin()) == '#') continue;

		temp_ss.str(""); temp_ss.clear(); temp_ss << temp_str;
		temp_ss >> temp_vert.posis.lon >> temp_vert.posis.lat >> temp_id;
		temp_vert.posis.rad = defaultR;
		temp_vert.set(temp_vert.posis);
		listVert.push_back(temp_vert);
		listIds.push_back(temp_id); //list id 可能不止一个 但是我们暂时只用一个就行了
	}
	infile.close();
	return 0;
}

int mshInterpolate::initRef(char* para)
{
	//首先匹配预定义类型
	if (!strcmp(para,"WGS84"))
	{
		refr = pole_radius;
		refR = equator_radius;
	}
	else if (!strcmp(para,"EarthRadius"))
	{
		refr = EarthRadius;
		refR = EarthRadius;
	}
	//匹配格式化参数
	else if (2 == sscanf(para,"%lf/%lf",&refr,&refR))
	{
		if (refr <= 0 || refR <= 0)
		{
			cout << BOLDRED << "Error ==> " << RESET << "fail to initial reference system" << endl;
			return -1;
		}
	}
	else
	{
		cout << BOLDRED << "Error ==> " << RESET << "fail to initial reference system" << endl;
		return -1;
	}
	return 0;
}

void mshInterpolate::interpolate_topo()
{
	vertex cross_point;
	vertex temp_vert[3];
	cpoint temp_nor;

	listVal.resize(listVert.size(),MAX_DBL);
	for (int i = 0; i < listVert.size(); i++)
	{
		for (int j = 0; j < 3; j++)
		{
			temp_vert[j] = mshVert[mshTri[listIds[i]].vec[j]];
			temp_vert[j].posis.rad += vertVal[mshTri[listIds[i]].vec[j]];
			temp_vert[j].set(temp_vert[j].posis);
		}

		temp_nor = cross(temp_vert[1].posic-temp_vert[0].posic,temp_vert[2].posic-temp_vert[0].posic);
		cross_point.posic = lineOnPlane(temp_vert[0].posic,temp_nor,listVert[i].posic);
		cross_point.set(cross_point.posic);
		listVal[i] = cross_point.posis.rad - REF_r(cross_point.posis.lat,refr,refR);
	}
	return;
}

void mshInterpolate::interpolate_linear()
{
	vertex cross_point;
	vertex temp_vert[3];
	cpoint temp_nor;

	listVal.resize(listVert.size(),MAX_DBL);
	for (int i = 0; i < listVert.size(); i++)
	{
		for (int j = 0; j < 3; j++)
		{
			temp_vert[j] = mshVert[mshTri[listIds[i]].vec[j]];
		}

		temp_nor = cross(temp_vert[1].posic-temp_vert[0].posic,temp_vert[2].posic-temp_vert[0].posic);
		cross_point.posic = lineOnPlane(temp_vert[0].posic,temp_nor,listVert[i].posic);
		cross_point.set(cross_point.posic);
		listVal[i] = triInterp_area(cross_point.posic,temp_vert[0].posic,temp_vert[1].posic,temp_vert[2].posic,
			vertVal[mshTri[listIds[i]].vec[0]],
			vertVal[mshTri[listIds[i]].vec[1]],
			vertVal[mshTri[listIds[i]].vec[2]]);
	}
	return;
}

int mshInterpolate::outList()
{
	for (int i = 0; i < listVert.size(); i++)
	{
		cout << listVert[i].posis.lon << " " << listVert[i].posis.lat << " " << setprecision(16) << listVal[i] << endl;
	}
	return 0;
}
#endif