Merge pull request #16 from gridapapps/Code_refactoring

Code refactoring

Author: oriolcg

.github/workflows/ci.yml

@@ -0,0 +1,97 @@
+name: CI
+on: [push, pull_request]
+jobs:
+  test:
+    name: Tests ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.5'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
+
+        
+        
+  drivers:
+    name: Drivers ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.5'
+        os:
+          - ubuntu-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - run: |
+          julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
+            using Pkg; Pkg.instantiate()'
+      - run: |
+          julia --color=yes --project=. --check-bounds=yes --depwarn=error -e '
+            (1,) .== 1; cd("test"); include("runtests.jl")'
+  docs:
+    name: Documentation
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: '1.5'
+      - run: |
+          julia --project=docs -e '
+            using Pkg
+            Pkg.develop(PackageSpec(path=pwd()))
+            Pkg.instantiate()'
+#      - run: |
+#          julia --project=docs -e '
+#            using Documenter: doctest
+#            using Gridap
+#            doctest(Gridap)'
+      - run: julia --project=docs docs/make.jl
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

models/elasticFlag.json

@@ -101,7 +101,7 @@ function execute(problem::Problem{:analytical};kwargs...)
     d < 1.0e-8
   end
   oldcell_to_coods = get_cell_coordinates(trian)
-  oldcell_to_is_in = collect1d(apply(is_in,oldcell_to_coods))
+  oldcell_to_is_in = collect1d(lazy_map(is_in,oldcell_to_coods))
   incell_to_cell = findall(oldcell_to_is_in)
   outcell_to_cell = findall(collect(Bool, .! oldcell_to_is_in))
   model_solid = DiscreteModel(model,incell_to_cell)
@@ -132,16 +132,15 @@ function execute(problem::Problem{:analytical};kwargs...)
   # Quadratures
   println("Defining quadratures")
   order = _get_kwarg(:order,kwargs,2)
-  order = _get_kwarg(:order,kwargs,2)
-  quads = get_FSI_quadratures(Tₕ,order)
+  dTₕ = get_FSI_measures(Tₕ,order)
 
   # Test FE Spaces
   println("Defining FE spaces")
   Y_ST, X_ST, Y_FSI, X_FSI = get_FE_spaces(strategy,coupling,models,order,bconds,constraint=:zeromean)
 
   # Stokes problem for initial solution
   println("Defining Stokes operator")
-  op_ST = get_Stokes_operator([X_ST,Y_ST],strategy,Tₕ[:Ωf],quads[:Ωf],μ_f,fv_ST_Ωf(0.0))
+  op_ST = get_Stokes_operator(X_ST,Y_ST,strategy,dTₕ[:Ωf],μ_f,fv_ST_Ωf(0.0))
 
   # Setup equation parameters
   mesh_params = Dict{Symbol,Any}(
@@ -172,7 +171,7 @@ function execute(problem::Problem{:analytical};kwargs...)
 
   # FSI problem
   println("Defining FSI operator")
-  op_FSI = get_FSI_operator([X_FSI,Y_FSI],coupling,strategy,Tₕ,quads,params)
+  op_FSI = get_FSI_operator(X_FSI,Y_FSI,coupling,strategy,Tₕ,dTₕ,params)
 
   # Setup output files
   folderName = "fsi-results"
@@ -184,53 +183,53 @@ function execute(problem::Problem{:analytical};kwargs...)
 
   # Solve Stokes problem
   @timeit "ST problem" begin
-  println("Defining Stokes solver")
-  xh = solve(op_ST)
-  if(is_vtk)
-    writePVD(filePath, Tₕ[:Ωf], [(xh, 0.0)])
+    println("Defining Stokes solver")
+    xh = solve(op_ST)
+    if(is_vtk)
+      writePVD(filePath, Tₕ[:Ωf], [(xh, 0.0)])
+    end
   end
-end
 
-# Compute Stokes solution L2-norm
-l2(w) = w⋅w
-eu_ST = u(0.0) - restrict(xh[1],Tₕ[:Ωf])
-ev_ST = v(0.0) - restrict(xh[2],Tₕ[:Ωf])
-ep_ST = p(0.0) - restrict(xh[3],Tₕ[:Ωf])
-eul2_ST = sqrt(sum( integrate(l2(eu_ST),Tₕ[:Ωf],quads[:Ωf]) ))
-evl2_ST = sqrt(sum( integrate(l2(ev_ST),Tₕ[:Ωf],quads[:Ωf]) ))
-epl2_ST = sqrt(sum( integrate(l2(ep_ST),Tₕ[:Ωf],quads[:Ωf]) ))
-println("Stokes L2-norm u: ", eul2_ST)
-println("Stokes L2-norm v: ", evl2_ST)
-println("Stokes L2-norm p: ", epl2_ST)
-@test eul2_ST < 1.0e-10
-@test evl2_ST < 1.0e-10
-@test epl2_ST < 1.0e-10
+  # Compute Stokes solution L2-norm
+  l2(w) = w⋅w
+  eu_ST = u(0.0) - xh[1]
+  ev_ST = v(0.0) - xh[2]
+  ep_ST = p(0.0) - xh[3]
+  eul2_ST = sqrt(∑( ∫(l2(eu_ST))dTₕ[:Ωf] ))
+  evl2_ST = sqrt(∑( ∫(l2(ev_ST))dTₕ[:Ωf] ))
+  epl2_ST = sqrt(∑( ∫(l2(ep_ST))dTₕ[:Ωf] ))
+  println("Stokes L2-norm u: ", eul2_ST)
+  println("Stokes L2-norm v: ", evl2_ST)
+  println("Stokes L2-norm p: ", epl2_ST)
+  @test eul2_ST < 1.0e-10
+  @test evl2_ST < 1.0e-10
+  @test epl2_ST < 1.0e-10
 
-# Solve FSI problem
-@timeit "FSI problem" begin
-println("Defining FSI solver")
-xh0 = interpolate_everywhere([u(0.0),v(0.0),p(0.0)],X_FSI(0.0))
-nls = NLSolver(
-show_trace = true,
-method = :newton,
-linesearch = BackTracking(),
-ftol = 1.0e-10,
-iterations = 50
-)
-odes =  ThetaMethod(nls, dt, 0.5)
-solver = TransientFESolver(odes)
-xht = solve(solver, op_FSI, xh0, t0, tf)
-end
+  # Solve FSI problem
+  @timeit "FSI problem" begin
+  println("Defining FSI solver")
+  xh0 = interpolate_everywhere([u(0.0),v(0.0),p(0.0)],X_FSI(0.0))
+  nls = NLSolver(
+  show_trace = true,
+  method = :newton,
+  linesearch = BackTracking(),
+  ftol = 1.0e-10,
+  iterations = 50
+  )
+  odes =  ThetaMethod(nls, dt, 0.5)
+  solver = TransientFESolver(odes)
+  xht = solve(solver, op_FSI, xh0, t0, tf)
+  end
 
-# Compute outputs
-out_params = Dict(
-:u=>u,
-:v=>v,
-:p=>p,
-:filePath=>filePath,
-:is_vtk=>is_vtk
-)
-output = computeOutputs(xht,Tₕ,quads,strategy,out_params)
+  # Compute outputs
+  out_params = Dict(
+  :u=>u,
+  :v=>v,
+  :p=>p,
+  :filePath=>filePath,
+  :is_vtk=>is_vtk
+  )
+  output = computeOutputs(xht,Tₕ,dTₕ,strategy,out_params)
 
 end
 
@@ -305,7 +304,7 @@ function get_FE_spaces(problem::Problem{:analytical},strategy::WeakForms.MeshStr
   )
 end
 
-function computeOutputs(xht,Tₕ,quads,strategy,params)
+function computeOutputs(xht,Tₕ,dTₕ,strategy,params)
 
   # Unpack parameters
   u = params[:u]
@@ -326,18 +325,18 @@ function computeOutputs(xht,Tₕ,quads,strategy,params)
       println("============================")
 
       # Compute errors
-      eu = u(t) - restrict(xh[1],Tₕ[:Ω])
-      ev = v(t) - restrict(xh[2],Tₕ[:Ω])
-      ep = p(t) - restrict(xh[3],Tₕ[:Ω])
-      eul2 = sqrt(sum( integrate(l2(eu),Tₕ[:Ω],quads[:Ω] )))
-      evl2 = sqrt(sum( integrate(l2(ev),Tₕ[:Ω],quads[:Ω] )))
-      epl2 = sqrt(sum( integrate(l2(ep),Tₕ[:Ω],quads[:Ω] )))
+      eu = u(t) - xh[1]
+      ev = v(t) - xh[2]
+      ep = p(t) - xh[3]
+      eul2 = sqrt(∑( ∫(l2(eu))dTₕ[:Ω] ))
+      evl2 = sqrt(∑( ∫(l2(ev))dTₕ[:Ω] ))
+      epl2 = sqrt(∑( ∫(l2(ep))dTₕ[:Ω] ))
 
       # Write to PVD
       if(is_vtk)
-        uh = restrict(xh[1],Tₕ[:Ω])
-        vh = restrict(xh[2],Tₕ[:Ω])
-        ph = restrict(xh[3],Tₕ[:Ω])
+        uh = xh[1]
+        vh = xh[2]
+        ph = xh[3]
         pvd[t] = createvtk(
         Tₕ[:Ω],
         filePath * "_$t.vtu",

models/elasticFlagFine.json

@@ -82,8 +82,8 @@ function execute(problem::Problem{:elasticFlag}; kwargs...)
   println("Defining discrete model")
   modelName = _get_kwarg(:model,kwargs,"../models/elasticFlag.json")
   model = DiscreteModelFromFile(modelName)
-  model_solid = DiscreteModel(model,"solid")
-  model_fluid = DiscreteModel(model,"fluid")
+  model_solid = DiscreteModel(model,tags="solid")
+  model_fluid = DiscreteModel(model,tags="fluid")
   models = Dict(:Ω => model, :Ωf => model_fluid, :Ωs => model_solid)
 
   # Triangulations
@@ -93,15 +93,15 @@ function execute(problem::Problem{:elasticFlag}; kwargs...)
   # Quadratures
   println("Defining quadratures")
   order = _get_kwarg(:order,kwargs,2)
-  quads = get_FSI_quadratures(Tₕ,order)
+  dTₕ = get_FSI_measures(Tₕ,order)
 
   # Test FE Spaces
   println("Defining FE spaces")
   Y_ST, X_ST, Y_FSI, X_FSI = get_FE_spaces(strategy,coupling,models,order,bconds)
 
   # Stokes problem for initial solution
   println("Defining Stokes operator")
-  op_ST = get_Stokes_operator([X_ST,Y_ST],strategy,Tₕ[:Ωf],quads[:Ωf],μ_f,f(0.0))
+  op_ST = get_Stokes_operator(X_ST,Y_ST,strategy,dTₕ[:Ωf],μ_f,f(0.0))
 
   # Setup equation parameters
   mesh_params = Dict{Symbol,Any}(
@@ -132,7 +132,7 @@ function execute(problem::Problem{:elasticFlag}; kwargs...)
 
   # FSI problem
   println("Defining FSI operator")
-  op_FSI = get_FSI_operator([X_FSI,Y_FSI],coupling,strategy,Tₕ,quads,params)
+  op_FSI = get_FSI_operator(X_FSI,Y_FSI,coupling,strategy,Tₕ,dTₕ,params)
 
   # Setup output files
   folderName = "fsi-results"
@@ -178,7 +178,7 @@ out_params = Dict{Symbol,Any}(
   :filePath=>filePath,
   :is_vtk=>is_vtk,
   )
-output = computeOutputs(xh0,xht,coupling,strategy,models,Tₕ,quads,out_params)
+output = computeOutputs(xh0,xht,coupling,strategy,models,Tₕ,dTₕ,out_params)
 
 end
 
@@ -248,12 +248,12 @@ function computeOutputs(xh0,xht,coupling::WeakForms.Coupling,strategy::WeakForms
   end
 
   ## Surface triangulation
-  trian_Γc = BoundaryTriangulation(models[:Ω], "cylinder")
-  quad_Γc = CellQuadrature(trian_Γc, bdegree)
-  n_Γc = get_normal_vector(trian_Γc)
+  Γc = BoundaryTriangulation(models[:Ω], tags="cylinder")
+  dΓc = Measure(Γc, bdegree)
+  n_Γc = get_normal_vector(Γc)
 
   # Aux function
-  traction_boundary(n,u,v,p) = n ⋅ WeakForms.Pf_dev(μ,u,v)  + WeakForms.Pf_vol(u,p) * n
+  traction_boundary(n,u,v,p) = n ⋅ WeakForms.Pᵥ_Ωf(μ,u,v)  + WeakForms.Pₚ_Ωf(u,p) * n
   function traction_closure(coupling,n,u,v,p)
     if typeof(coupling) == WeakForms.Coupling{:weak}
       traction_boundary(n,u,v,p).⁺
@@ -285,30 +285,13 @@ function computeOutputs(xh0,xht,coupling::WeakForms.Coupling,strategy::WeakForms
     uh = xh[uvpindex[1]]
     vh = xh[uvpindex[2]]
     ph = xh[uvpindex[3]]
-    uh_Γc = restrict(uh, trian_Γc)
-    vh_Γc = restrict(vh, trian_Γc)
-    ph_Γc = restrict(ph, trian_Γc)
-    uhn_Γc = restrict(uhn, trian_Γc)
-    vhn_Γc = restrict(vhn, trian_Γc)
-    phn_Γc = restrict(phn, trian_Γc)
-    uhθ_Γc = θ*uh_Γc + (1.0-θ)*uhn_Γc
-    vhθ_Γc = θ*vh_Γc + (1.0-θ)*vhn_Γc
-    phθ_Γc = θ*ph_Γc + (1.0-θ)*phn_Γc
-    uh_Γi = restrict(uh, Tₕ[:Γi])
-    vh_Γi = restrict(vh, Tₕ[:Γi])
-    ph_Γi = restrict(ph, Tₕ[:Γi])
-    uhn_Γi = restrict(uhn, Tₕ[:Γi])
-    vhn_Γi = restrict(vhn, Tₕ[:Γi])
-    phn_Γi = restrict(phn, Tₕ[:Γi])
-    uhθ_Γi = θ*uh_Γi + (1.0-θ)*uhn_Γi
-    vhθ_Γi = θ*vh_Γi + (1.0-θ)*vhn_Γi
-    phθ_Γi = θ*ph_Γi + (1.0-θ)*phn_Γi
+    xθ(xₙ₊₁,xₙ) = θ*xₙ₊₁ + (1-θ)*xₙ
 
     # Integrate on the cylinder
-    FDc, FLc = sum(integrate(traction_boundary(n_Γc,uhθ_Γc,vhθ_Γc,phθ_Γc), trian_Γc, quad_Γc))
+    FDc, FLc = ∑(∫(traction_boundary(n_Γc,xθ(uh,uhn),xθ(vh,vhn),xθ(ph,phn)) )dΓc )
 
     # Integrate on the interface
-    FDi, FLi = sum(integrate(traction_interface(n_Γi,uhθ_Γi,vhθ_Γi,phθ_Γi), Tₕ[:Γi], quads[:Γi]))
+    FDi, FLi = ∑(∫(traction_boundary(n_Γc,xθ(uh,uhn),xθ(vh,vhn),xθ(ph,phn)) )dΓc )
 
     FD = FDc + FDi
     FL = FLc + FLi