"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Functions/F_Raytracing.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

F_Raytracing.cpp  (getdp-3.4.0-source.tgz):F_Raytracing.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include <math.h> #include <math.h>
#include "GetDPConfig.h" #include "GetDPConfig.h"
#include "ProData.h" #include "ProData.h"
#include "F.h" #include "F.h"
#include "Message.h" #include "Message.h"
skipping to change at line 29 skipping to change at line 29
void F_DiamondPhase(F_ARG) void F_DiamondPhase(F_ARG)
{ {
Message::Error("F_DiamondPhase requires the GSL"); Message::Error("F_DiamondPhase requires the GSL");
} }
#else #else
#include <gsl/gsl_vector.h> #include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h> #include <gsl/gsl_multiroots.h>
struct f_context struct f_context {
{
double x1, y1; double x1, y1;
}; };
static int f(const gsl_vector *ts, void *param, gsl_vector *f) static int f(const gsl_vector *ts, void *param, gsl_vector *f)
{ {
struct f_context * c = (struct f_context *)param; struct f_context *c = (struct f_context *)param;
double t = gsl_vector_get(ts, 0); double t = gsl_vector_get(ts, 0);
double s = gsl_vector_get(ts, 1); double s = gsl_vector_get(ts, 1);
double x = c->x1; double x = c->x1;
double y = c->y1; double y = c->y1;
gsl_vector_set(f, 0, cos(t) - s*cos(2*t) - x); gsl_vector_set(f, 0, cos(t) - s * cos(2 * t) - x);
gsl_vector_set(f, 1, sin(t) - s*sin(2*t) - y); gsl_vector_set(f, 1, sin(t) - s * sin(2 * t) - y);
return GSL_SUCCESS; return GSL_SUCCESS;
} }
static int df(const gsl_vector *ts, void* param, gsl_matrix *j) static int df(const gsl_vector *ts, void *param, gsl_matrix *j)
{ {
double t = gsl_vector_get(ts, 0); double t = gsl_vector_get(ts, 0);
double s = gsl_vector_get(ts, 1); double s = gsl_vector_get(ts, 1);
double j1dt = -sin(t) + s*2*sin(2*t); double j1dt = -sin(t) + s * 2 * sin(2 * t);
double j2dt = cos(t) - 2*s*cos(2*t); double j2dt = cos(t) - 2 * s * cos(2 * t);
double j1ds = -cos(2*t); double j1ds = -cos(2 * t);
double j2ds = -sin(2*t); double j2ds = -sin(2 * t);
gsl_matrix_set(j, 0, 0, j1dt); gsl_matrix_set(j, 0, 0, j1dt);
gsl_matrix_set(j, 1, 1, j2ds); gsl_matrix_set(j, 1, 1, j2ds);
gsl_matrix_set(j, 1, 0, j2dt); gsl_matrix_set(j, 1, 0, j2dt);
gsl_matrix_set(j, 0, 1, j1ds); gsl_matrix_set(j, 0, 1, j1ds);
return GSL_SUCCESS; return GSL_SUCCESS;
} }
static int fdf(const gsl_vector *uv, void *param, gsl_vector *func, gsl_matrix * static int fdf(const gsl_vector *uv, void *param, gsl_vector *func,
jac) gsl_matrix *jac)
{ {
f(uv, param, func); f(uv, param, func);
df(uv, param, jac); df(uv, param, jac);
return GSL_SUCCESS; return GSL_SUCCESS;
} }
static int newton(gsl_multiroot_function_fdf FDF, double *u, double *v) static int newton(gsl_multiroot_function_fdf FDF, double *u, double *v)
{ {
const int MAX_ITER = 25; const int MAX_ITER = 25;
const gsl_multiroot_fdfsolver_type* TYPE = gsl_multiroot_fdfsolver_gnewton; const gsl_multiroot_fdfsolver_type *TYPE = gsl_multiroot_fdfsolver_gnewton;
int iter = 0, status; int iter = 0, status;
gsl_multiroot_fdfsolver* solver = gsl_multiroot_fdfsolver_alloc(TYPE, 2); gsl_multiroot_fdfsolver *solver = gsl_multiroot_fdfsolver_alloc(TYPE, 2);
/* u, v contains initial guess */ /* u, v contains initial guess */
gsl_vector *X = gsl_vector_alloc(2); gsl_vector *X = gsl_vector_alloc(2);
gsl_vector_set(X, 0, *u); gsl_vector_set(X, 0, *u);
gsl_vector_set(X, 1, *v); gsl_vector_set(X, 1, *v);
gsl_multiroot_fdfsolver_set(solver, &FDF, X); gsl_multiroot_fdfsolver_set(solver, &FDF, X);
do { do {
iter++; iter++;
status = gsl_multiroot_fdfsolver_iterate(solver); status = gsl_multiroot_fdfsolver_iterate(solver);
*u = gsl_vector_get(solver->x, 0); *u = gsl_vector_get(solver->x, 0);
*v = gsl_vector_get(solver->x, 1); *v = gsl_vector_get(solver->x, 1);
if(*v < 0 || *v > 15 || fabs(*u) > 7){ if(*v < 0 || *v > 15 || fabs(*u) > 7) {
status= GSL_FAILURE; status = GSL_FAILURE;
break; break;
} }
status = gsl_multiroot_test_residual(solver->f, 1.e-12); status = gsl_multiroot_test_residual(solver->f, 1.e-12);
} while(status == GSL_CONTINUE && iter < MAX_ITER); } while(status == GSL_CONTINUE && iter < MAX_ITER);
gsl_multiroot_fdfsolver_free(solver); gsl_multiroot_fdfsolver_free(solver);
gsl_vector_free(X); gsl_vector_free(X);
if(status == GSL_SUCCESS) if(status == GSL_SUCCESS)
skipping to change at line 126 skipping to change at line 126
double x = A->Val[0], y = A->Val[1]; double x = A->Val[0], y = A->Val[1];
struct f_context context = {x, y}; struct f_context context = {x, y};
gsl_multiroot_function_fdf FDF; gsl_multiroot_function_fdf FDF;
if(x > 0 && y < 1 && y > -1) { if(x > 0 && y < 1 && y > -1) {
V->Val[0] = x; V->Val[0] = x;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
if(x > 0){ if(x > 0) {
tau[1] = sqrt(x * x + y * y); tau[1] = sqrt(x * x + y * y);
if(y > 0){ if(y > 0) { initGuess = (atan2(y, -x) + 3.14 / 2) / 2; }
initGuess = (atan2(y, -x) + 3.14 / 2) / 2; else {
}
else{
initGuess = (atan2(y, -x) - 3.14 / 2) / 2; initGuess = (atan2(y, -x) - 3.14 / 2) / 2;
} }
} }
else{ else {
tau[1] = sqrt(x * x + y * y) - 1; tau[1] = sqrt(x * x + y * y) - 1;
initGuess = atan2(y, x); initGuess = atan2(y, x);
} }
if(fabs(x) < 1 && fabs(y) > 6.5){ if(fabs(x) < 1 && fabs(y) > 6.5) {
if(y < 0){ if(y < 0) { initGuess = initGuess - 3.14 / 8; }
initGuess = initGuess - 3.14 / 8; else {
}
else{
initGuess = initGuess + 3.14 / 8; initGuess = initGuess + 3.14 / 8;
} }
} }
tau[0] = initGuess; tau[0] = initGuess;
if(tau[1] == 0){ if(tau[1] == 0) {
V->Val[0] = x; V->Val[0] = x;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
FDF.f = &f; FDF.f = &f;
FDF.df = &df; FDF.df = &df;
FDF.fdf = &fdf; FDF.fdf = &fdf;
FDF.n = 2; FDF.n = 2;
FDF.params = &context; FDF.params = &context;
if(!newton(FDF, &tau[0], &tau[1])) if(!newton(FDF, &tau[0], &tau[1]))
Message::Error("Newton did not converge: %lf, %lf \n", tau[0], tau[1]); Message::Error("Newton did not converge: %lf, %lf \n", tau[0], tau[1]);
/* now we just go on to calculate the phase from this */ /* now we just go on to calculate the phase from this */
phase = cos(tau[0]) + tau[1]; phase = cos(tau[0]) + tau[1];
if(phase > abs(13)){ if(phase > abs(13)) { phase = 13; }
phase = 13;
}
V->Val[0] = phase; V->Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
} }
void F_DiamondPhase(F_ARG) void F_DiamondPhase(F_ARG)
{ {
double x, y, phase, theta, xtrans, ytrans; double x, y, phase, theta, xtrans, ytrans;
x = A->Val[0]; x = A->Val[0];
skipping to change at line 197 skipping to change at line 191
phase = -x; phase = -x;
V-Val[0] = phase; V-Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
*/ */
x = -x; /* just a temp investigation */ x = -x; /* just a temp investigation */
/*partition up the space into a couple of pieces*/ /*partition up the space into a couple of pieces*/
if(x >= 0 && (y-.1 <= 1 && y+.1 >= -1)){ if(x >= 0 && (y - .1 <= 1 && y + .1 >= -1)) {
V->Val[0] = x; V->Val[0] = x;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
/* /*
if( x <= 0 && (y>=-1 && y<=0) ) if( x <= 0 && (y>=-1 && y<=0) )
{ {
phase = -y-1 + (-x+(-y-1)); phase = -y-1 + (-x+(-y-1));
V->Val[0] = phase; V->Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
*/ */
/*check to see if the point is in the cone made by the x-corner*/ /*check to see if the point is in the cone made by the x-corner*/
xtrans = x + 1; xtrans = x + 1;
ytrans = y; ytrans = y;
theta = atan2(ytrans,xtrans); theta = atan2(ytrans, xtrans);
if( theta >= 3.14/2 || theta <= -3.14/2){ if(theta >= 3.14 / 2 || theta <= -3.14 / 2) {
phase = -1 + sqrt( pow(xtrans,2.0) + pow(ytrans,2.0) ); phase = -1 + sqrt(pow(xtrans, 2.0) + pow(ytrans, 2.0));
V->Val[0] = phase; V->Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
/*check to see if the point is in the upper corner cone*/ /*check to see if the point is in the upper corner cone*/
xtrans = x; xtrans = x;
ytrans = y - 1; ytrans = y - 1;
theta = atan2(ytrans,xtrans); theta = atan2(ytrans, xtrans);
if( theta <= 3.14/2 && theta >= 0 ){ if(theta <= 3.14 / 2 && theta >= 0) {
phase = sqrt( pow(x,2.0) + pow(ytrans,2.0) ); phase = sqrt(pow(x, 2.0) + pow(ytrans, 2.0));
V->Val[0] = phase; V->Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
/*lower corner cone*/ /*lower corner cone*/
xtrans = x; xtrans = x;
ytrans = y + 1; ytrans = y + 1;
theta = atan2(ytrans,xtrans); theta = atan2(ytrans, xtrans);
if( theta >= -3.14/2 && theta <= 0 ){ if(theta >= -3.14 / 2 && theta <= 0) {
phase = sqrt( pow(x,2.0) + pow(ytrans,2.0) ); phase = sqrt(pow(x, 2.0) + pow(ytrans, 2.0));
V->Val[0] = phase; V->Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
return; return;
} }
/*the point must be in one of the two reflections caused by the /*the point must be in one of the two reflections caused by the
sides facing the incoming wave*/ sides facing the incoming wave*/
/* xtrans = x; /* xtrans = x;
ytrans = y; */ ytrans = y; */
if(y<0){ if(y < 0) {
/* xtrans = (x+y+1)/2; /* xtrans = (x+y+1)/2;
ytrans = (xtrans)-1; ytrans = (xtrans)-1;
phase = xtrans + sqrt( pow(x-xtrans,2.0) + pow(y-ytrans,2.0) ); */ phase = xtrans + sqrt( pow(x-xtrans,2.0) + pow(y-ytrans,2.0) ); */
phase = -x + ( -y + (x-1) ); phase = -x + (-y + (x - 1));
} }
else{ else {
/* xtrans = (1+x-y)/2; /* xtrans = (1+x-y)/2;
ytrans = xtrans-1; ytrans = xtrans-1;
phase = xtrans + sqrt( pow(x-xtrans,2.0) + pow(-y-ytrans,2.0) ); */ phase = xtrans + sqrt( pow(x-xtrans,2.0) + pow(-y-ytrans,2.0) ); */
phase = -x + ( y - (1-x) ); phase = -x + (y - (1 - x));
} }
V->Val[0] = phase; V->Val[0] = phase;
V->Type = SCALAR; V->Type = SCALAR;
} }
#endif #endif
 End of changes. 24 change blocks. 
50 lines changed or deleted 43 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)