tetra_vz/main_msh.cpp

#include "tetra_vz.h"
#include "string"
#include "sstream"
#include "iostream"
#include "fstream"

std::stringstream str2ss(std::string in)
{
	std::stringstream tmp_ss;
	tmp_ss.clear();
	tmp_ss.str(in);
	return tmp_ss;
}

int main(int argc, char *argv[])
{
	int pnum, tnum; // Size of the vertexes and tetrahedrons
	point3d *modelp = nullptr; // vertex class
	tetrahedron *modelt = nullptr; // tetrahedron class

	// file object
	std::ifstream infile("cube.msh");

	int tmp_i;
	int tmp_tnum;
	std::string line_str;
	std::stringstream line_ss;
	while (getline(infile, line_str))
	{
		if (line_str == "$Nodes") // Read vertexes
		{
			// Get vertex size firstly
			getline(infile, line_str);
			line_ss = str2ss(line_str);
			line_ss >> pnum;
			// Allocate memories
			modelp = new point3d [pnum];
			for (int i = 0; i < pnum; i++)
			{
				getline(infile, line_str);
				line_ss = str2ss(line_str);
				line_ss >> tmp_i >> modelp[i].x >> modelp[i].y >> modelp[i].z;
			}
		}
		else if (line_str == "$Elements") // Read tetrahedral indexes
		{
			// Get tetrahedron's sizes
			getline(infile, line_str);
			line_ss = str2ss(line_str);
			line_ss >> tmp_tnum;
			tnum = tmp_tnum;

			for (int i = 0; i < tmp_tnum; i++)
			{
				getline(infile, line_str);
				line_ss = str2ss(line_str);
				line_ss >> tmp_i >> tmp_i;
				if (tmp_i != 4) // The tag for tetrahedron is 4
				{
					tnum--;
				}
				else break;
			}

			// Allocate memories
			modelt = new tetrahedron [tnum];

			line_ss = str2ss(line_str);

			for (int j = 0; j < 5; j++) // Skip the first 5 numbers
				line_ss >> tmp_i;
			for (int j = 0; j < 4; j++)
			{
				line_ss >> tmp_i;
				modelt[0].vec_ptr[j] = &modelp[tmp_i-1];
			}

			for (int i = 1; i < tnum; i++)
			{
				getline(infile, line_str);
				line_ss = str2ss(line_str);

				for (int j = 0; j < 5; j++) // Skip the first 5 numbers
					line_ss >> tmp_i;
				for (int j = 0; j < 4; j++)
				{
					line_ss >> tmp_i;
					modelt[i].vec_ptr[j] = &modelp[tmp_i-1];
				}
			}
		}
	}

	infile.close();

	// Output for hunman validation
	std::clog << "node number = " << pnum << std::endl;
	std::clog << "tetra number = " << tnum << std::endl;

	// Set tetrahedron's density and initialize it
	for (int i = 0; i < tnum; i++)
	{
		modelt[i].rho = 1.0;
		modelt[i].initialize_tensors();
	}

	// Declare observation stations
	double xmin = 0, dx = 2.5; // 41
	double ymin = 0, dy = 2.5; // 41
	point3d obs[1681]; // 41*41
	for (int i = 0; i < 41; i++)
	{
		for (int j = 0; j < 41; j++)
		{
			obs[i*41+j].x = xmin + dx*j;
			obs[i*41+j].y = ymin + dy*i;
			obs[i*41+j].z = 0.0;
		}
	}

	// Calculate the gravity value
	double sum;
	for (int i = 0; i < 1681; i++)
	{
		std::cout << obs[i].x << " " << obs[i].y << " ";

		sum = 0.0;
		for (int t = 0; t < tnum; t++)
		{
			sum += tetra_vz(&modelt[t], &obs[i]);
		}
		std::cout << sum << std::endl;
	}

	if (modelp != nullptr) delete[] modelp;
	if (modelt != nullptr) delete[] modelt;
	return 0;
}