Non-linear multi field problems working

src/FESpaces/FEOperators.jl

@@ -32,7 +32,7 @@ import Gridap.FESpaces: TrialFESpace
 import Gridap.FESpaces: TestFESpace
 import Gridap.LinearSolvers: LinearSolver
 
-# We use duck typing in this module
+# We use duck typing in this module for the following types:
 
 """
 Arguments annotated with this type have to implement the following queries:

src/MultiField/MultiFEFunctions.jl

@@ -14,6 +14,7 @@ import Base: getindex
 import Base: iterate
 import Base: zero
 import Gridap.FESpaces: free_dofs
+import Gridap.FESpaces: FEFunction
 
 struct MultiFEFunction
   fields::Vector{<:FEFunction}
@@ -28,6 +29,11 @@ function MultiFEFunction(
   MultiFEFunction(fields,free_dofs_all_fields)
 end
 
+function FEFunction(
+  fespaces::MultiFESpace, free_dofs_all_fields::AbstractVector)
+  MultiFEFunction(free_dofs_all_fields,fespaces)
+end
+
 free_dofs(self::MultiFEFunction) = self.free_dofs_all_fields
 
 length(self::MultiFEFunction) = length(self.fields)

src/MultiField/MultiFEOperators.jl

@@ -7,6 +7,7 @@ using Gridap.Geometry
 using Gridap.CellQuadratures
 using Gridap.FESpaces
 using Gridap.FEOperators: LinearFEOperator
+using Gridap.FEOperators: NonLinearFEOperator
 using ..MultiFESpaces
 using ..MultiAssemblers
 using ..MultiFEFunctions
@@ -25,4 +26,17 @@ function LinearFEOperator(
   LinearFEOperator(biform,liform,V,U,assem,trian,quad)
 end
 
+function NonLinearFEOperator(
+  res::Function,
+  jac::Function,
+  testfesp::Vector{<:FESpaceWithDirichletData},
+  trialfesp::Vector{<:FESpaceWithDirichletData},
+  assem::MultiAssembler,
+  trian::Triangulation{Z},
+  quad::CellQuadrature{Z}) where Z
+  V = MultiFESpace(testfesp)
+  U = MultiFESpace(trialfesp)
+  NonLinearFEOperator(res,jac,V,U,assem,trian,quad)
+end
+
 end # module MultiFEOperators

test/MultiFieldTests/MultiFEFunctionsTests.jl

@@ -51,4 +51,8 @@ zh = zero(U)
 @test isa(zh,MultiFEFunction)
 @test free_dofs(zh) == zeros(num_free_dofs(U))
 
+zh = FEFunction(U,free_dofs(zh))
+@test isa(zh,MultiFEFunction)
+@test free_dofs(zh) == zeros(num_free_dofs(U))
+
 end # module MultiFEFunctionsTests

test/MultiFieldTests/MultiNonLinearFEOperatorsTests.jl

@@ -0,0 +1,113 @@
+include("../../src/MultiField/MultiFEOperators.jl")
+
+module MultiFEOperatorsTests
+
+using Gridap.FESpaces
+using Gridap.Assemblers
+using Gridap.FEOperators
+using Gridap.LinearSolvers
+using Gridap.NonLinearSolvers
+
+using Test
+using Gridap
+using Gridap.Geometry
+using Gridap.CellMaps
+using Gridap.Geometry.Cartesian
+using Gridap.FieldValues
+using Gridap.CellQuadratures
+using Gridap.CellIntegration
+using Gridap.Vtkio
+
+import Gridap: gradient
+
+# Define manufactured functions
+u1fun(x) = x[1] + x[2]
+u2fun(x) = x[1] - x[2]
+
+u1fun_grad(x) = VectorValue(1.0,1.0)
+u2fun_grad(x) = VectorValue(1.0,-1.0)
+
+gradient(::typeof(u1fun)) = u1fun_grad
+gradient(::typeof(u2fun)) = u2fun_grad
+
+b1fun(x) = u2fun(x) -(3.0*x[1]+x[2]+1.0)
+b2fun(x) = 0.0
+
+νfun(x,u1) = (u1+1.0)*x[1]
+
+# Construct the discrete model
+model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))
+
+# Construct the FEspace
+order = 1
+diritag = "boundary"
+fespace = ConformingFESpace(Float64,model,order,diritag)
+
+# Define test and trial
+V1 = TestFESpace(fespace)
+V2 = V1
+V = [V1, V2]
+
+U1 = TrialFESpace(fespace,u1fun)
+U2 = TrialFESpace(fespace,u2fun)
+U = [U1, U2]
+
+# Define integration mesh and quadrature
+trian = Triangulation(model)
+quad = CellQuadrature(trian,order=2)
+
+# Define cell field describing the source term
+b1field = CellField(trian,b1fun)
+b2field = CellField(trian,b2fun)
+
+# Define a solution dependent material parameter
+ν(u1) = CellField(trian,νfun,u1)
+
+# Define residual and jacobian
+a(u,v,du) = inner(∇(v[1]),ν(u[1])*∇(du[1])) + inner(v[1],du[2]) + inner(∇(v[2]),∇(du[2]))
+b(v) = inner(v[1],b1field) + inner(v[2],b2field)
+
+res(u,v) = a(u,v,u) - b(v)
+jac(u,v,du) = a(u,v,du) # + inner(∇(v[1]),ν(du[1])*∇(u[1]))
+
+# Define Assembler
+assem = SparseMatrixAssembler(V,U)
+
+# Define the FEOperator
+op = NonLinearFEOperator(res,jac,V,U,assem,trian,quad)
+
+# Define the FESolver
+ls = LUSolver()
+tol = 1.e-10
+maxiters = 20
+nls = NewtonRaphsonSolver(ls,tol,maxiters)
+solver = NonLinearFESolver(nls)
+
+# Solve!
+uh = solve(solver,op)
+
+# Define exact solution and error
+u1 = CellField(trian,u1fun)
+e1 = u1 - uh[1]
+
+u2 = CellField(trian,u2fun)
+e2 = u2 - uh[2]
+
+# Define norms to measure the error
+l2(u) = inner(u,u)
+h1(u) = inner(∇(u),∇(u)) + l2(u)
+
+# Compute errors
+e1l2 = sqrt(sum( integrate(l2(e1),trian,quad) ))
+e1h1 = sqrt(sum( integrate(h1(e1),trian,quad) ))
+
+e2l2 = sqrt(sum( integrate(l2(e2),trian,quad) ))
+e2h1 = sqrt(sum( integrate(h1(e2),trian,quad) ))
+
+@test e1l2 < 1.e-8
+@test e1h1 < 1.e-8
+
+@test e2l2 < 1.e-8
+@test e2h1 < 1.e-8
+
+end

test/MultiFieldTests/runtests.jl

@@ -9,5 +9,6 @@ using Test
 @testset "MultiFEFunctions" begin include("MultiFEFunctionsTests.jl") end
 @testset "MultiFEBases" begin include("MultiFEBasesTests.jl") end
 @testset "MultiFEOperators" begin include("MultiFEOperatorsTests.jl") end
+@testset "MultiNonLinearFEOperators" begin include("MultiNonLinearFEOperatorsTests.jl") end
 
 end # module MultiFieldTests