diff --git a/lib/potential/gkernel_triangle2d.cpp b/lib/potential/gkernel_triangle2d.cpp index b26d516..fca6b8e 100644 --- a/lib/potential/gkernel_triangle2d.cpp +++ b/lib/potential/gkernel_triangle2d.cpp @@ -42,6 +42,8 @@ typedef void (*gkernel_triangle2d_ptr)(gctl::matrix &out_kernel, const g void gkernel_triangle2d_vz(gctl::matrix &out_kernel, const gctl::array &ele, const gctl::array &obsp, gctl::verbose_type_e verbose); +void gkernel_triangle2d_vx(gctl::matrix &out_kernel, const gctl::array &ele, + const gctl::array &obsp, gctl::verbose_type_e verbose); void gkernel_triangle2d_vzx(gctl::matrix &out_kernel, const gctl::array &ele, const gctl::array &obsp, gctl::verbose_type_e verbose); void gkernel_triangle2d_vzz(gctl::matrix &out_kernel, const gctl::array &ele, @@ -65,6 +67,9 @@ void gctl::gkernel(matrix &out_kernel, const array &ele, case Vz: triangle_kernel = gkernel_triangle2d_vz; break; + case Vx: + triangle_kernel = gkernel_triangle2d_vx; + break; case Tzx: triangle_kernel = gkernel_triangle2d_vzx; break; @@ -92,6 +97,8 @@ typedef void (*gobser_triangle2d_ptr)(gctl::array &out_obs, const gctl:: void gobser_triangle2d_vz(gctl::array &out_obs, const gctl::array &ele, const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose); +void gobser_triangle2d_vx(gctl::array &out_obs, const gctl::array &ele, + const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose); void gobser_triangle2d_vzx(gctl::array &out_obs, const gctl::array &ele, const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose); void gobser_triangle2d_vzz(gctl::array &out_obs, const gctl::array &ele, @@ -116,6 +123,9 @@ void gctl::gobser(array &out_obs, const array &ele, const ar case Vz: triangle_obser = gobser_triangle2d_vz; break; + case Vx: + triangle_obser = gobser_triangle2d_vx; + break; case Tzx: triangle_obser = gobser_triangle2d_vzx; break; @@ -130,6 +140,7 @@ void gctl::gobser(array &out_obs, const array &ele, const ar } double gkernel_triangle2d_vz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr); +double gkernel_triangle2d_vx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr); double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr); double gkernel_triangle2d_vzz_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr); @@ -174,7 +185,31 @@ void gkernel_triangle2d_vz(gctl::matrix &out_kernel, const gctl::array &out_kernel, const gctl::array &ele, + const gctl::array &obsp, gctl::verbose_type_e verbose) +{ + int i, j; + int o_size = obsp.size(); + int e_size = ele.size(); + + out_kernel.resize(o_size, e_size); + + gctl::progress_bar bar(o_size, "gkernel_vx"); + for (i = 0; i < o_size; i++) + { + if (verbose == gctl::FullMsg) bar.progressed(i); + else if (verbose == gctl::ShortMsg) bar.progressed_simple(i); + +#pragma omp parallel for private (j) schedule(guided) + for (j = 0; j < e_size; j++) + { + out_kernel[i][j] = gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i)); } } return; @@ -198,7 +233,7 @@ void gkernel_triangle2d_vzx(gctl::matrix &out_kernel, const gctl::array< #pragma omp parallel for private (j) schedule(guided) for (j = 0; j < e_size; j++) { - out_kernel.at(i, j) = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)); + out_kernel[i][j] = gkernel_triangle2d_vzx_sig(ele.get(j), obsp.get(i)); } } return; @@ -222,7 +257,7 @@ void gkernel_triangle2d_vzz(gctl::matrix &out_kernel, const gctl::array< #pragma omp parallel for private (j) schedule(guided) for (j = 0; j < e_size; j++) { - out_kernel.at(i, j) = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)); + out_kernel[i][j] = gkernel_triangle2d_vzz_sig(ele.get(j), obsp.get(i)); } } return; @@ -246,7 +281,31 @@ void gobser_triangle2d_vz(gctl::array &out_obs, const gctl::array &out_obs, const gctl::array &ele, + const gctl::array &obsp, const gctl::array &rho, gctl::verbose_type_e verbose) +{ + int i, j; + int o_size = obsp.size(); + int e_size = ele.size(); + + out_obs.resize(o_size, 0.0); + + gctl::progress_bar bar(e_size, "gobser_vx"); + for (j = 0; j < e_size; j++) + { + if (verbose == gctl::FullMsg) bar.progressed(j); + else if (verbose == gctl::ShortMsg) bar.progressed_simple(j); + +#pragma omp parallel for private (i) schedule(guided) + for (i = 0; i < o_size; i++) + { + out_obs[i] += gkernel_triangle2d_vx_sig(ele.get(j), obsp.get(i)) * rho[j]; } } return; @@ -270,7 +329,7 @@ void gobser_triangle2d_vzx(gctl::array &out_obs, const gctl::array &out_obs, const gctl::arrayvert[n]->y - op_ptr->y)*(ele_ptr->vert[(n+1)%3]->x - op_ptr->x) - + (ele_ptr->vert[(n+1)%3]->y - op_ptr->y)*(ele_ptr->vert[n]->x - op_ptr->x); + Ab = pow(ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y,2) + + pow(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x,2); + A = Aa/Ab; + + Ba = atan2(ele_ptr->vert[n]->y - op_ptr->y, ele_ptr->vert[n]->x - op_ptr->x); + Bb = atan2(ele_ptr->vert[(n+1)%3]->y - op_ptr->y, ele_ptr->vert[(n+1)%3]->x - op_ptr->x); + B = (ele_ptr->vert[(n+1)%3]->y - ele_ptr->vert[n]->y)*(Ba - Bb); + + Ca = pow(ele_ptr->vert[(n+1)%3]->x - op_ptr->x,2) + pow(ele_ptr->vert[(n+1)%3]->y - op_ptr->y,2); + Cb = pow(ele_ptr->vert[n]->x - op_ptr->x,2) + pow(ele_ptr->vert[n]->y - op_ptr->y,2); + C = 0.5*(ele_ptr->vert[(n+1)%3]->x - ele_ptr->vert[n]->x)*(log(Ca) - log(Cb)); + + sum += A*(B+C); + } + return -2.0e+8*GCTL_G0*sum; +} + double gkernel_triangle2d_vzx_sig(gctl::triangle2d *ele_ptr, gctl::point2dc *op_ptr) { double sum;