fmm-examples/fmm_2d_triangle/fmm_2d_triangle.cpp

#include "iostream"
#include "fstream"
#include "sstream"
#include "string.h"
#include "cmath"
#include "vector"
#include "iomanip"
#include "algorithm"
#include "ctime"

#include "gctl/core.h"
#include "gctl/geometry.h"
#include "gctl/io.h"

using namespace std;

struct vertex;
struct ele;

struct source : public gctl::point2dc
{
	double rad;
	double slow;
};

struct vertex : public gctl::vertex2dc
{
	int tag = 0; //0 = far away, 1 = close, 2 = active
	double slow; // slow are constant within a element. a mean value must be set in local update
	double time = 1e+30;
	double syn_time = 1e+30; //synthetic direct arrive time (only for error test)
	vector<vertex*> v_neigh; //neighbor vertex unknown amount
	vector<ele*> e_neigh; //neighbor element unknown amount
	vector<int> ve_order; // 此顶点在相接三角形中的顶点序号(0、1或者2)，与e_neigh中的三角形一一对应，应该能提高计算效率
};

struct ele
{
	int id, tag = 0; //0 = not calculated, 1 = already calculated
	vertex* vec_ptr[3];
	double slow;
	double area;
	double edge_l[3];
};

double local_update_triangle(double area, double length_ab, double length_ac, 
	double length_bc, double a_time, double b_time, double slow)
{
	///< 沿边ac或bc传播的走时值
	double direct_time = GCTL_MIN(a_time + slow*length_ac, 
		b_time + slow*length_bc);

	double u = b_time - a_time;
	double w = slow*slow*length_ab*length_ab - u*u;
	if (w < 0) return direct_time;

	w = sqrt(w);
	double rho0 = 2.0*area/length_ab;
	double ksi0 = sqrt(length_ac*length_ac - rho0*rho0)/length_ab;
	double ksi  = ksi0 - u*rho0/(length_ab*w);

	if (ksi >= 1 || ksi <= 0) return direct_time;
	else
	{
		double trial_time = a_time + u*ksi0 + w*rho0/length_ab;
		if (trial_time > a_time && trial_time > b_time)
			return GCTL_MIN(trial_time, direct_time);
		else return direct_time;
	}
}

void local_update(vertex* vert_ptr)
{
	int j;
	double time1, time2;
	// loop over all neighboring triangles. If any one have two known vertice, update time
	for (int i = 0; i < vert_ptr->e_neigh.size(); i++)
	{
		if (vert_ptr->e_neigh[i]->tag == 0)
		{
			j = vert_ptr->ve_order[i];
			if (vert_ptr->e_neigh[i]->vec_ptr[(j+1)%3]->tag == 2 &&
				vert_ptr->e_neigh[i]->vec_ptr[(j+2)%3]->tag == 2)
			{
				time1 = vert_ptr->e_neigh[i]->vec_ptr[(j+1)%3]->time;
				time2 = vert_ptr->e_neigh[i]->vec_ptr[(j+2)%3]->time;

				if (time1 <= time2)
					vert_ptr->time = GCTL_MIN(vert_ptr->time, local_update_triangle(
						vert_ptr->e_neigh[i]->area, vert_ptr->e_neigh[i]->edge_l[j],
						vert_ptr->e_neigh[i]->edge_l[(j+2)%3], vert_ptr->e_neigh[i]->edge_l[(j+1)%3],
						time1, time2, vert_ptr->e_neigh[i]->slow));
				else
					vert_ptr->time = GCTL_MIN(vert_ptr->time, local_update_triangle(
						vert_ptr->e_neigh[i]->area, vert_ptr->e_neigh[i]->edge_l[j],
						vert_ptr->e_neigh[i]->edge_l[(j+1)%3], vert_ptr->e_neigh[i]->edge_l[(j+2)%3],
						time2, time1, vert_ptr->e_neigh[i]->slow));
				vert_ptr->e_neigh[i]->tag = 1;
			}
		}
	}
	return;
}

class FMM_2D_TRIANGLE
{
public:
	FMM_2D_TRIANGLE(){}
	~FMM_2D_TRIANGLE(){}
	// read mesh files generated by triangle
	int ReadFiles(char* filename);
	// output a Gmsh (.msh) file
	int OutMsh(char* filename, bool tag);
	// initialize node slowness. slowness of the source will be used if default slowness is zero.
	void InitNodeSlowness(double default_slow = 0.0);
	// add abnormal slowness
	void SetRectangeSlowness(double ab_slow,double xmin,double xmax,double ymin,double ymax);
	// initialize triangle's slowness
	void InitEleSlowness();
	// calculate synthetic direct arrive time and fmm time
	void CalculateSolution();
	// set source parameters
	int set_init_source(double x,double y,double r,double s);
private:
	int node_num_, ele_num_;
	vector<vertex> nodes_;
	vector<ele> elements_;

	vector<vertex*> nodes_ptr_;
	source init_source_;
};

int FMM_2D_TRIANGLE::ReadFiles(char* filename)
{
	ifstream infile;
	string temp_str;
	stringstream temp_ss;

	char full_name[1024];
	//read node file first
	strcpy(full_name,filename);
	strcat(full_name,".node");
	infile.open(full_name);
	if (!infile)
	{
		cerr << "file open error: " << full_name << endl;
		return -1;
	}
	else
	{
		clog << "reading node file:\t" << full_name;

		getline(infile,temp_str);
		gctl::str2ss(temp_str, temp_ss);
		temp_ss >> node_num_;

		nodes_.resize(node_num_);
		for (int i = 0; i < node_num_; i++)
		{
			getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> nodes_[i].id >> nodes_[i].x >> nodes_[i].y;
		}
		infile.close();

		clog << " ... done!" << endl;
	}

	//read element file
	int tmp_ids[3];
	strcpy(full_name,filename);
	strcat(full_name,".ele");
	infile.open(full_name);
	if (!infile)
	{
		cerr << "file open error: " << full_name << endl;
		return -1;
	}
	else
	{
		clog << "reading element file:\t" << full_name;

		getline(infile,temp_str);
		gctl::str2ss(temp_str, temp_ss);
		temp_ss >> ele_num_;

		elements_.resize(ele_num_);
		for (int i = 0; i < ele_num_; i++)
		{
			getline(infile,temp_str);
			gctl::str2ss(temp_str, temp_ss);
			temp_ss >> elements_[i].id >> tmp_ids[0] >> tmp_ids[1] >> tmp_ids[2];
			elements_[i].vec_ptr[0] = &nodes_[tmp_ids[0]];
			elements_[i].vec_ptr[1] = &nodes_[tmp_ids[1]];
			elements_[i].vec_ptr[2] = &nodes_[tmp_ids[2]];
		}
		infile.close();

		clog << " ... done!" << endl;
	}

	// reorder the vertice sequence to anti-clockwise
	gctl::point2dc v01, v02;
	vertex *tmp_ptr;
	for (int i = 0; i < ele_num_; i++)
	{
		v01.x = elements_[i].vec_ptr[1]->x - elements_[i].vec_ptr[0]->x;
		v01.y = elements_[i].vec_ptr[1]->y - elements_[i].vec_ptr[0]->y;

		v02.x = elements_[i].vec_ptr[2]->x - elements_[i].vec_ptr[0]->x;
		v02.y = elements_[i].vec_ptr[2]->y - elements_[i].vec_ptr[0]->y;

		if (v01.x*v02.y - v01.y*v02.x < 0)
		{
			tmp_ptr = elements_[i].vec_ptr[1];
			elements_[i].vec_ptr[1] = elements_[i].vec_ptr[2];
			elements_[i].vec_ptr[2] = tmp_ptr;
		}
	}

	//sort neighboring vertice and triangles for a vertex
	double p;
	for (int i = 0; i < ele_num_; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			// opposite edge length of vertex j
			elements_[i].edge_l[j] = gctl::distance(*elements_[i].vec_ptr[(j+1)%3], *elements_[i].vec_ptr[(j+2)%3]);
		}

		// calculate triangle's area
		p = 0.5*(elements_[i].edge_l[0] + elements_[i].edge_l[1] + elements_[i].edge_l[2]);
		elements_[i].area = sqrt(p*(p-elements_[i].edge_l[0])*(p-elements_[i].edge_l[1])*(p-elements_[i].edge_l[2]));

		for (int j = 0; j < 3; j++)
		{
			// add current element to each vertex element's neighboring list
			elements_[i].vec_ptr[j]->e_neigh.push_back(&elements_[i]);
			elements_[i].vec_ptr[j]->ve_order.push_back(j);
			// add the other two vertice to current vertex's node neighboring list
			// this operation will create duplicate enters for nodes's neighbors
			elements_[i].vec_ptr[j]->v_neigh.push_back(elements_[i].vec_ptr[(j+1)%3]);
			elements_[i].vec_ptr[j]->v_neigh.push_back(elements_[i].vec_ptr[(j+2)%3]);
		}
	}

	// remove duplicate enters for neighboring vertice
	vector<vertex*>::iterator pos;
	for (int i = 0; i < node_num_; i++)
	{
		sort(nodes_[i].v_neigh.begin(),nodes_[i].v_neigh.end()); //对顶点序列由小到大排序
		pos = unique(nodes_[i].v_neigh.begin(),nodes_[i].v_neigh.end()); //获取重复序列开始的位置
		nodes_[i].v_neigh.erase(pos,nodes_[i].v_neigh.end()); //删除重复点
	}

	//output statistics to screen
	clog << "node number:\t" << node_num_ << endl;
	clog << "element number:\t" << ele_num_ << endl;

	return 0;
}

int FMM_2D_TRIANGLE::OutMsh(char* filename, bool tag)
{
	ofstream outfile;
	if (tag)
		strcat(filename,"_linear.msh");
	else
		strcat(filename,".msh");

	outfile.open(filename);
	if (!outfile)
	{
		cerr << "file open error: " << filename << endl;
		return -1;
	}
	else clog << "writing .msh file:\t" << filename << endl;

	outfile<<"$MeshFormat"<<endl<<"2.2 0 8"<<endl<<"$EndMeshFormat"<<endl<<"$Nodes"<<endl<<node_num_<<endl;
	// set the first vertex index to 1 to avoid a display bug in Gmsh
	for (int i = 0; i < node_num_; i++)
	{
		outfile << nodes_[i].id + 1 << " " << setprecision(16) << nodes_[i].x << " " << nodes_[i].y << " 0" << endl;
	}
	outfile<<"$EndNodes"<<endl;
	outfile<<"$Elements"<<endl<<ele_num_<<endl;
	for (int i = 0; i < ele_num_; i++)
	{
		outfile << elements_[i].id + 1 <<" 2 1 0";
		for (int j = 0; j < 3; j++)
			outfile << " " << elements_[i].vec_ptr[j]->id + 1;
		outfile << endl;
	}
	outfile << "$EndElements"<< endl;
	if (tag)
	{
		int tmp_node_num = 0;
		for (int i = 0; i < node_num_; i++)
		{
			if (nodes_[i].time < 1e+10)
				tmp_node_num++;
		}
		outfile <<"$NodeData"<< endl;
		outfile<<1<<endl<<"\"fmm arrive time (ms)\""<<endl<<1<<endl<<0.0<<endl
		<<3<<endl<<0<<endl<<1<<endl<<tmp_node_num<<endl;
		for (int i = 0; i < node_num_; i++)
		{
			if (nodes_[i].time < 1e+10)
			{
				outfile << nodes_[i].id + 1 << " " << setprecision(16) << nodes_[i].time << endl;
			}
		}
		outfile<<"$EndNodeData"<< endl;

		outfile <<"$NodeData"<< endl;
		outfile<<1<<endl<<"\"synthetic time (ms)\""<<endl<<1<<endl<<0.0<<endl
		<<3<<endl<<0<<endl<<1<<endl<<tmp_node_num<<endl;
		for (int i = 0; i < node_num_; i++)
		{
			if (nodes_[i].time < 1e+10) // this condition seems unnecessary, we will keep it here until we know better
			{
				outfile << nodes_[i].id + 1 << " " << setprecision(16) << nodes_[i].syn_time << endl;
			}
		}
		outfile<<"$EndNodeData"<< endl;

		// 这里我们顺带统计一下误差的平均值和标准差
		vector<double> tmp_err(tmp_node_num, 0.0);
		double mean_err = 0.0;
		outfile <<"$NodeData"<< endl;
		outfile<<1<<endl<<"\"fmm - synthetic time (ms)\""<<endl<<1<<endl<<0.0<<endl
		<<3<<endl<<0<<endl<<1<<endl<<tmp_node_num<<endl;
		for (int i = 0; i < node_num_; i++)
		{
			if (nodes_[i].time < 1e+10)
			{
				tmp_err[i] = nodes_[i].time - nodes_[i].syn_time;
				mean_err += tmp_err[i];

				outfile << nodes_[i].id + 1 << " " << setprecision(16) << tmp_err[i] << endl;
			}
		}
		outfile<<"$EndNodeData"<< endl;

		mean_err /= 1.0*tmp_node_num;
		double rms_err = 0.0;
		for (int i = 0; i < tmp_node_num; i++)
		{
			rms_err += pow(tmp_err[i] - mean_err,2);
		}
		clog << "mean error  = " << mean_err << endl;
		clog << "rms of error= " << sqrt(rms_err/tmp_node_num) << endl;

		outfile <<"$ElementData"<< endl;
		outfile<<1<<endl<<"\"model slowness (s/m)\""<<endl<<1<<endl<<0.0<<endl
		<<3<<endl<<0<<endl<<1<<endl<<ele_num_<<endl;
		for (int i = 0; i < ele_num_; i++)
		{
			outfile << elements_[i].id + 1 << " " << setprecision(16) << elements_[i].slow << endl;
		}
		outfile<<"$EndElementData"<< endl;
	}
	outfile.close();
	return 0;
}

void FMM_2D_TRIANGLE::InitNodeSlowness(double default_slow)
{
	// if no background slowness is given. 
	// set all node's slowness as the same as source's slowness.
	if (default_slow == 0.0)
	{
		for (int i = 0; i < node_num_; i++)
		{
			nodes_[i].slow = init_source_.slow;
		}
	}
	else
	{
		for (int i = 0; i < node_num_; i++)
		{
			nodes_[i].slow = default_slow;
		}
	}
	return;
}

// set slowness within a rectangular area.
void FMM_2D_TRIANGLE::SetRectangeSlowness(double ab_slow,double xmin,double xmax,double ymin,double ymax)
{
	for (int i = 0; i < node_num_; i++)
	{
		if (nodes_[i].x >= xmin && nodes_[i].x <= xmax &&
			nodes_[i].y >= ymin && nodes_[i].y <= ymax)
		{
			nodes_[i].slow = ab_slow;
		}
	}
	return;
}

// set element's slowness as the mean value of three vertice's slowness.
void FMM_2D_TRIANGLE::InitEleSlowness()
{
	for (int i = 0; i < ele_num_; i++)
	{
		elements_[i].slow = (elements_[i].vec_ptr[0]->slow + elements_[i].vec_ptr[1]->slow + elements_[i].vec_ptr[2]->slow)/3.0;
	}
	return;
}

void FMM_2D_TRIANGLE::CalculateSolution()
{
	time_t start_time, stop_time;
	start_time = time(NULL);
	clog << "calculating first arrive time ... ";

	//calculate synthetic direct arrival time for all nodes
	//this calculation is only valid for testing when there is no abnormal slowness or obstacles.
	for (int i = 0; i < node_num_; i++)
	{
		nodes_[i].syn_time = gctl::distance(nodes_[i], init_source_) * init_source_.slow;
	}

	//initialize source nodes and close nodes
	//set a vertex's tag as 'active' if it is inside the initial radius
	//and calculate node's time as direct arrival time
	for (int i = 0; i < node_num_; i++)
	{
		if (gctl::distance(nodes_[i], init_source_) <= init_source_.rad)
		{
			nodes_[i].tag = 2;
			nodes_[i].time = gctl::distance(nodes_[i], init_source_) * init_source_.slow;
		}
	}

	//set all close vertice's tag to 'close' and add them to wave front list
	for (int i = 0; i < node_num_; i++)
	{
		if (nodes_[i].tag == 2)
		{
			for (int j = 0; j < nodes_[i].v_neigh.size(); j++)
			{
				if (nodes_[i].v_neigh[j]->tag == 0)
				{
					nodes_[i].v_neigh[j]->tag = 1;
					nodes_ptr_.push_back(nodes_[i].v_neigh[j]);
				}
			}
		}
	}

	//assign initial tags for elements
	//set an element as calculated only if all its three vectice are activated.
	for (int i = 0; i < ele_num_; i++)
	{
		if (elements_[i].vec_ptr[0]->tag == 2 &&
			elements_[i].vec_ptr[1]->tag == 2 &&
			elements_[i].vec_ptr[2]->tag == 2)
		{
			elements_[i].tag = 1;
		}
	}

	//calculate trial time for all close nodes
	for (int i = 0; i < nodes_ptr_.size(); i++)
	{
		local_update(nodes_ptr_[i]);
	}

	gctl::heap_sort<vertex*> time_sorter;
	//marching forward and updating the close nodes set
	while (!nodes_ptr_.empty())
	{
		// heap sort close nodes pointers to put a node at the first place if it has the smallest time
		std::sort(nodes_ptr_.begin(), nodes_ptr_.end(), [](vertex *a, vertex *b)->bool{
			if (a->time < b->time) return true; return false;
		});

		//change the first node's tag to 2 and update it's neighbor's time if it is not active (tag != 2)
		//this step is needed as a vertex may possess smaller time once it's neighbor time is updated.
		nodes_ptr_[0]->tag = 2;
		for (int i = 0; i < nodes_ptr_[0]->v_neigh.size(); i++)
		{
			// if a neighboring node is not added to the wave front list.
			// Add it to the list and update it's time
			if (nodes_ptr_[0]->v_neigh[i]->tag == 0)
			{
				nodes_ptr_[0]->v_neigh[i]->tag = 1;
				local_update(nodes_ptr_[0]->v_neigh[i]);
				nodes_ptr_.push_back(nodes_ptr_[0]->v_neigh[i]);
			}
			// already in the wave front list, update it's time
			else if (nodes_ptr_[0]->v_neigh[i]->tag == 1)
			{
				local_update(nodes_ptr_[0]->v_neigh[i]);
			}
		}
		// remove the first node from the wave front list
		nodes_ptr_.erase(nodes_ptr_.begin());
	}

	stop_time = time(NULL);
	clog << stop_time - start_time << " s" << endl;
	return;
}

int FMM_2D_TRIANGLE::set_init_source(double x, double y, double r, double s)
{
	if (s <=0 || r <= 0)
	{
		cerr << "source initialization error!" << endl;
		return -1;
	}
	else
	{
		init_source_.x = x; init_source_.y = y;
		init_source_.slow = s; init_source_.rad = r;
		return 0;
	}
}

/*******************************************main function here*****************************/
// argv[1] -> running mode
// argv[2] -> name of the input triangle file
// argv[3] -> source parameter
// argv[4] -> abnormal slownesses separated by commas

int main(int argc, char* argv[])
{
	FMM_2D_TRIANGLE instance;
	if (!strcmp(argv[1],"convert"))
	{
		if(!instance.ReadFiles(argv[2]))
			instance.OutMsh(argv[2],false);
	}
	else if (!strcmp(argv[1],"calculate"))
	{
		double s_x, s_y, s_rad, s_slow;
		if (4 != sscanf(argv[3], "%lf/%lf/%lf/%lf", &s_x, &s_y, &s_rad, &s_slow))
		{
			cerr << "wrong source parameter." << endl;
			return 0;
		}

		if(instance.set_init_source(s_x, s_y, s_rad, s_slow)) return 0;
		if(instance.ReadFiles(argv[2])) return 0;
		//set node slowness here
		instance.InitNodeSlowness();

		//add abnormal slowness
		double ab_slow, ab_xmin, ab_xmax, ab_ymin, ab_ymax;
		string all_ab_para, old_ab_para;
		string tmp_str;
		stringstream tmp_ss;
		if (argc >= 5)
		{
			old_ab_para = argv[4];
			gctl::replace_all(all_ab_para, old_ab_para, ",", " ");
			gctl::str2ss(all_ab_para, tmp_ss);
			while (tmp_ss >> tmp_str)
			{
				if (5 != sscanf(tmp_str.c_str(), "%lf/%lf/%lf/%lf/%lf", 
					&ab_slow, &ab_xmin, &ab_xmax, &ab_ymin, &ab_ymax))
				{
					cerr << "wrong source parameter." << endl;
					return 0;
				}
				instance.SetRectangeSlowness(ab_slow, ab_xmin, ab_xmax, ab_ymin, ab_ymax);
			}
		}

		// set element slowness here
		instance.InitEleSlowness();
		instance.CalculateSolution();
		if(instance.OutMsh(argv[2],true)) return 0;
	}
	return 0;
}