/** ![Victor Gustave Robin could have looked like Peter Gustav Lejeune Dirichlet (or Claude Louis Marie Henri Navier)](https://alchetron.com/cdn/peter-gustav-lejeune-dirichlet-4db41bd3-b41c-4a34-b3b9-80d9a762eff-resize-750.jpeg) # Mixed boundary conditions From [Yong Hui's page](../yonghui/smalltest/robin.c) we know that for a boundary condition of the form, $$a \mathtt{s} + b \frac{\partial{\mathtt{s}}}{\partial \mathbf{n}} = c,$$ The discretized version for $2b\neq -a\Delta$ reads, $$s [ \mathtt{ghost} ] = \left[ \frac{2 c \Delta}{2b + a\Delta} \right]+ \Bigg \langle \frac{2b - a\Delta }{2b + a\Delta} \Bigg \rangle s[\quad]$$ The ansatz is that this may be expressed as a *linear mix* of a dirichlet and a Neumann condition. $$s[\mathtt{ghost}] = A\cdot \mathtt{dirichlet}(D) + B\cdot\mathtt{neumann}(N).$$ Expanding the macros gives: $$s[\mathtt{ghost}] = 2AD - As[\quad]+ B\Delta N + Bs[ \quad]$$ $$= \left[2AD + B\Delta N\right]+ \langle B-A\rangle s[\quad].$$ Each term at the rhs forms a seperate equation. Such that we have four unkowns ($A,B,D,N$) and only two equations. Such underdetermined system may have many solutions. Lets check if a solution exist for $A = 1$ and $N = 0$. The equation for the square-bracketed term then becomes: $$2D = \frac{2 c \Delta}{2b + a\Delta},\rightarrow$$ $$D = \frac{c \Delta}{2b + a\Delta}.$$ Next, for the second, angle-bracketed-term equation, $$B - 1 = \frac{2b - a\Delta }{2b + a\Delta} \rightarrow$$ $$B = \frac{2b - a\Delta }{2b + a\Delta} + 1.$$ Wrapping it up (in many braces): */ #define robin(a,b,c) ((dirichlet ((c)*Delta/(2*(b) + (a)*Delta))) + ((neumann (0))*((2*(b) - (a)*Delta)/(2*(b) + (a)*Delta) + 1.))) /** Mind the devide-by-zero risk. ## A test [Vyaas Gururajan](https://groups.google.com/g/basilisk-fr/c/-YSthke6VBo) coded a neat test case. Special care must be taken, because it is not always critical. ~~~gnuplot Looks OK set xlabel 'x' set ylabel 'u' set size square set grid plot 'out' u 1:3 w l lw 6 t 'Analytical Solution', '' w l lw 3 t 'Approximate Solution' ~~~ */ #include "grid/multigrid1D.h" #include "poisson.h" /** Here is the analaytical solution taken from [Professor Fitzpatrick's notes](http://farside.ph.utexas.edu/teaching/329/lectures/node66.html): */ double alphaL = 2., betaL = -1., gammaL = -3.; double alphaR = 1., betaR = 3., gammaR = -2.; double analyticalU (double x) { double d = alphaL*alphaR + alphaL*betaR - betaL*alphaR; double g = (gammaL*(alphaR + betaR) - betaL*(gammaR - (alphaR+betaR)/3.))/d; double h = (alphaL*(gammaR - (alphaR + betaR)/3.) - gammaL*alphaR)/d; return (g + h*x + sq(x)/2. - sq(sq(x))/6.); } scalar u[], rhs[]; u[left] = robin (alphaL, -betaL, gammaL); //Mind the sign! u[right] = robin (alphaR, betaR, gammaR); int main() { init_grid (1 << 7); foreach() rhs[] = 1. - 2.*sq(x); poisson (u, rhs, tolerance = 1e-9); foreach() printf("%g %g %g\n",x, u[], analyticalU(x)); }