using Gridap
domain = (0,1,0,1)
partition = (100,100)
model = CartesianDiscreteModel(domain,partition)

order = 1

V = TestFESpace(
  reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64},
  conformity=:HDiv, model=model, dirichlet_tags=[5,6])

Q = TestFESpace(
  reffe=:QLagrangian, order=order, valuetype=Float64,
  conformity=:L2, model=model)

uD = VectorValue(0.0,0.0)
U = TrialFESpace(V,uD)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

neumanntags = [8,]
btrian = BoundaryTriangulation(model,neumanntags)
bquad = CellQuadrature(btrian,degree)

const kinv1 = TensorValue(1.0,0.0,0.0,1.0)
const kinv2 = TensorValue(100.0,90.0,90.0,100.0)
@law function σ(x,u)
   if ((abs(x[1]-0.5) <= 0.1) && (abs(x[2]-0.5) <= 0.1))
      return kinv2*u
   else
      return kinv1*u
   end
end

px = get_physical_coordinate(trian)

function a(x,y)
   v, q = y
   u, p = x
   v*σ(px,u) - (∇*v)*p + q*(∇*u)
end

nb = get_normal_vector(btrian)
h = -1.0
function b_ΓN(y)
  v, q = y
  (v*nb)*h
end

t_Ω = LinearFETerm(a,trian,quad)
t_ΓN = FESource(b_ΓN,btrian,bquad)
op = AffineFEOperator(X,Y,t_Ω,t_ΓN)
xh = solve(op)
uh, ph = xh

writevtk(trian,"darcyresults",cellfields=["uh"=>uh,"ph"=>ph])