tmp update

fmm_2d_triangle/fmm_2d_triangle.cpp

@@ -17,16 +17,9 @@ 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
@@ -111,23 +104,20 @@ public:
 	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();
+	void InitEleSlowness(double in_slow);
 	// calculate synthetic direct arrive time and fmm time
-	void CalculateSolution();
+	void CalculateSolution(double bkg_slow);
 	// set source parameters
-	int set_init_source(double x,double y,double r,double s);
+	void set_init_source(const gctl::point2dc &loc);
 private:
-	int node_num_, ele_num_;
+	size_t node_num_, ele_num_, src_id_;
 	vector<vertex> nodes_;
 	vector<ele> elements_;
 
 	vector<vertex*> nodes_ptr_;
-	source init_source_;
 };
 
 int FMM_2D_TRIANGLE::ReadFiles(char* filename)
@@ -364,52 +354,32 @@ int FMM_2D_TRIANGLE::OutMsh(char* filename, bool tag)
 	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++)
+	gctl::point2dc cen;
+	for (int i = 0; i < ele_num_; i++)
 	{
-		if (nodes_[i].x >= xmin && nodes_[i].x <= xmax &&
-			nodes_[i].y >= ymin && nodes_[i].y <= ymax)
+		cen = 1.0/3.0 * (*elements_[i].vec_ptr[0] + *elements_[i].vec_ptr[1] + *elements_[i].vec_ptr[2]);
+		if (cen.x >= xmin && cen.x <= xmax && cen.y >= ymin && cen.y <= ymax)
 		{
-			nodes_[i].slow = ab_slow;
+			elements_[i].slow = ab_slow;
 		}
 	}
 	return;
 }
 
 // set element's slowness as the mean value of three vertice's slowness.
-void FMM_2D_TRIANGLE::InitEleSlowness()
+void FMM_2D_TRIANGLE::InitEleSlowness(double in_slow)
 {
 	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;
+		elements_[i].slow = in_slow;
 	}
 	return;
 }
 
-void FMM_2D_TRIANGLE::CalculateSolution()
+void FMM_2D_TRIANGLE::CalculateSolution(double bkg_slow)
 {
 	time_t start_time, stop_time;
 	start_time = time(NULL);
@@ -419,19 +389,7 @@ void FMM_2D_TRIANGLE::CalculateSolution()
 	//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;
-		}
+		nodes_[i].syn_time = gctl::distance(nodes_[i], nodes_[src_id_]) * bkg_slow;
 	}
 
 	//set all close vertice's tag to 'close' and add them to wave front list
@@ -450,18 +408,6 @@ void FMM_2D_TRIANGLE::CalculateSolution()
 		}
 	}
 
-	//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++)
 	{
@@ -504,19 +450,38 @@ void FMM_2D_TRIANGLE::CalculateSolution()
 	return;
 }
 
-int FMM_2D_TRIANGLE::set_init_source(double x, double y, double r, double s)
+void FMM_2D_TRIANGLE::set_init_source(const gctl::point2dc &loc)
 {
-	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;
-	}
+	// Find the closest node and make it the source
+    double dist, mini_dist = 1e+30;
+    for (size_t i = 0; i < node_num_; i++)
+    {
+        dist = gctl::distance(loc, nodes_[i]);
+        if (dist < mini_dist)
+        {
+            mini_dist = dist;
+            src_id_ = i;
+        }
+    }
+
+    // initiate the source time
+    nodes_[src_id_].time = 0.0;
+    nodes_[src_id_].tag = 2;
+    for (size_t i = 0; i < nodes_[src_id_].e_neigh.size(); i++)
+    {
+        nodes_[src_id_].e_neigh[i]->tag = 1;
+
+        for (size_t j = 0; j < 3; j++)
+        {
+            if (nodes_[src_id_].e_neigh[i]->vec_ptr[j]->id != src_id_)
+            {
+				dist = gctl::distance(nodes_[src_id_], *nodes_[src_id_].e_neigh[i]->vec_ptr[j]);
+                nodes_[src_id_].e_neigh[i]->vec_ptr[j]->time = nodes_[src_id_].e_neigh[i]->slow * dist;
+                nodes_[src_id_].e_neigh[i]->vec_ptr[j]->tag = 2;
+            }
+        }
+    }
+    return;
 }
 
 /*******************************************main function here*****************************/
@@ -535,17 +500,17 @@ int main(int argc, char* argv[])
 	}
 	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))
+		if(instance.ReadFiles(argv[2])) return 0;
+
+		double s_x, s_y, s_slow;
+		if (3 != sscanf(argv[3], "%lf/%lf/%lf", &s_x, &s_y, &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();
+		//set slowness here
+		instance.InitEleSlowness(s_slow);
+		instance.set_init_source(gctl::point2dc(s_x, s_y));
 
 		//add abnormal slowness
 		double ab_slow, ab_xmin, ab_xmax, ab_ymin, ab_ymax;
@@ -569,9 +534,7 @@ int main(int argc, char* argv[])
 			}
 		}
 
-		// set element slowness here
-		instance.InitEleSlowness();
-		instance.CalculateSolution();
+		instance.CalculateSolution(s_slow);
 		if(instance.OutMsh(argv[2],true)) return 0;
 	}
 	return 0;