// file nonsymelliptic.edp
//
// solution of $-a\Delta u + b\partial_1 u + cu = f$
// together with $u = u_0$ on $\Gamma$
//
// $a$, $b$, $c$ are positive constants
// $f = c x^2 + 2b x - c y^2$
// $u_0 = x^2 - y^2$

//
// NUMERICAL VALUES TO PROVIDE
//
int C0=100, C1=99, C2=98; // the labels
real a=.1, b=5., bt2 = 2.*b, c=10.; // the values on the boundaries
func f0=(c*x+bt2)*x - c*y*y; // the function f
func u00=x*x - y*y; // the function u_0

//
// DESCRIPTION OF THE BOUNDARY AND THE INTERFACE
//
// outer boundary ($\Gamma_0$):
border C00(t=0,2*pi){x=5*cos(t); y=5*sin(t); label=C0;}

// inner boundary ($\Gamma_1$):
border C11(t=0,1){ x=1+t;  y=3;      label=C1;}
border C12(t=0,1){ x=2;    y=3-6*t;  label=C1;}
border C13(t=0,1){ x=2-t;  y=-3;     label=C1;}
border C14(t=0,1){ x=1;    y=-3+6*t; label=C1;}

//
// FIGURES (I)
//
// boundary and interface:
plot (                C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20),
                    wait=true, ps="nonsymb.eps");

//
// MESH
//
mesh Th=buildmesh(    C00(50)
                    + C11(5)+C12(20)+C13(5)+C14(20));

//
// FIGURES (II)
//
// mesh:
plot(Th, wait=true, ps="nonsymTh.eps");

//
// THE FINITE ELEMENT SPACE AND FUNCTIONS
//
fespace Vh(Th,P1); Vh u,v;

//
// THE PROBLEM (FORMULATION AND SOLUTION)
//
problem nonsym(u,v) = int2d(Th)(a*(dx(u)*dx(v)+dy(u)*dy(v))
        + b*dx(u)*v + c*u*v)
        - int2d(Th)(f0*v)
        + on(C0,u=u00) + on(C1,u=u00);
nonsym;

//
// FIGURES (III)
//
// the solution:
plot(u, wait=true, value=true, fill=true, ps="nonsymsol.eps");