Laplacian

NEWS.md

@@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.15.1] - Unreleased
+
+### Added
+- Added support for Hessian and Laplacian operators. Only implemented for Finite Elements with an `AffineMap`.
+- Added a test for the solution of the Biharmonic equation.
+
 ## [0.15.0] - 2020-12-14
 
 This version is a major (backwards-incompatible) refactoring of the project which is not summarized here for the sake of brevity. Most of the functionality of v0.14.0 is available in v0.15.0, but possibly with a significantly different API. See [here](https://github.com/gridap/Tutorials/compare/v0.14.0...v0.15.0) the changes in the sources of the Gridap Tutorials between versions 0.14.0 and 0.15.0 to effectively see the major changes in the API.

src/CellData/CellFields.jl

@@ -274,6 +274,9 @@ function ∇∇(a::CellField)
   GenericCellField(h,get_triangulation(a),DomainStyle(a))
 end
 
+# This function has to be removed when ∇⋅∇(a) is implemented
+laplacian(a::CellField) = tr(∇∇(a))
+
 # Operations between CellField
 
 function evaluate!(cache,k::Operation,a::CellField...)

src/Fields/AffineMaps.jl

@@ -40,3 +40,32 @@ function gradient(h::AffineMap)
   ConstantField(h.gradient)
 end
 
+function push_∇∇(∇∇a::Field,ϕ::AffineMap)
+  # Assuming ϕ is affine map
+  Jt = ∇(ϕ)
+  Jt_inv = inv(Jt)
+  Operation(push_∇∇)(∇∇a, Jt_inv)
+end
+
+function push_∇∇(∇∇a::Number,Jt_inv::MultiValue{Tuple{D,D}} where D)
+  Jt_inv⋅Jt_inv⋅∇∇a
+end
+
+function lazy_map(
+  k::Broadcasting{typeof(push_∇∇)},
+  cell_∇∇a::AbstractArray,
+  cell_map::AbstractArray{<:AffineMap})
+  cell_Jt = lazy_map(∇,cell_map)
+  cell_invJt = lazy_map(Operation(inv),cell_Jt)
+  lazy_map(Broadcasting(Operation(push_∇∇)),cell_∇∇a,cell_invJt)
+end
+
+function lazy_map(
+  k::Broadcasting{typeof(push_∇∇)},
+  cell_∇∇at::LazyArray{<:Fill{typeof(transpose)}},
+  cell_map::AbstractArray{<:AffineMap})
+  cell_∇∇a = cell_∇∇at.args[1]
+  cell_∇∇b = lazy_map(k,cell_∇∇a,cell_map)
+  cell_∇∇bt = lazy_map(transpose,cell_∇∇b)
+  cell_∇∇bt
+end

src/Fields/ApplyOptimizations.jl

@@ -247,15 +247,19 @@ function lazy_map(k::Broadcasting{typeof(push_∇)},cell_∇a::AbstractArray,cel
   lazy_map(Broadcasting(Operation(⋅)),cell_invJt,cell_∇a)
 end
 
-function lazy_map(
-  k::Broadcasting{typeof(push_∇)},
-  cell_∇at::LazyArray{<:Fill{typeof(transpose)}},
-  cell_map::AbstractArray)
-
-  cell_∇a = cell_∇at.args[1]
-  cell_∇b = lazy_map(k,cell_∇a,cell_map)
-  cell_∇bt = lazy_map(transpose,cell_∇b)
-  cell_∇bt
+for op in (:push_∇,:push_∇∇)
+  @eval begin
+    function lazy_map(
+      k::Broadcasting{typeof($op)},
+      cell_∇at::LazyArray{<:Fill{typeof(transpose)}},
+      cell_map::AbstractArray)
+
+      cell_∇a = cell_∇at.args[1]
+      cell_∇b = lazy_map(k,cell_∇a,cell_map)
+      cell_∇bt = lazy_map(transpose,cell_∇b)
+      cell_∇bt
+    end
+  end
 end
 
 # Composing by the identity

src/Fields/DiffOperators.jl

@@ -62,8 +62,9 @@ const Δ = laplacian
     laplacian(f)
 """
 function laplacian(f)
-  g = gradient(f)
-  divergence(g)
+  # g = gradient(f)
+  # divergence(g)
+  tr(∇∇(f))
 end
 
 """
@@ -110,9 +111,8 @@ outer(f::Function,::typeof(∇)) = transpose(gradient(f))
     cross(∇,f)
 
 Equivalent to
-    
+
     curl(f)
 """
 cross(::typeof(∇),f::Field) = curl(f)
 cross(::typeof(∇),f::Function) = curl(f)
-

src/Fields/FieldArrays.jl

@@ -179,9 +179,13 @@ for T in (:(Point),:(AbstractVector{<:Point}))
   end
 end
 
-function gradient(a::LinearCombinationField)
-  fields = Broadcasting(∇)(a.fields)
-  LinearCombinationField(a.values,fields,a.column)
+for op in (:∇,:∇∇)
+  @eval begin
+    function $op(a::LinearCombinationField)
+      fields = Broadcasting($op)(a.fields)
+      LinearCombinationField(a.values,fields,a.column)
+    end
+  end
 end
 
 function linear_combination(a::AbstractMatrix{<:Number},b::AbstractVector{<:Field})
@@ -233,9 +237,13 @@ for T in (:(Point),:(AbstractVector{<:Point}))
   end
 end
 
-function evaluate!(cache,k::Broadcasting{typeof(∇)},a::LinearCombinationFieldVector)
-  fields = Broadcasting(∇)(a.fields)
-  LinearCombinationFieldVector(a.values,fields)
+for op in (:∇,:∇∇)
+  @eval begin
+    function evaluate!(cache,k::Broadcasting{typeof($op)},a::LinearCombinationFieldVector)
+      fields = Broadcasting($op)(a.fields)
+      LinearCombinationFieldVector(a.values,fields)
+    end
+  end
 end
 
 function get_children(n::TreeNode, a::LinearCombinationFieldVector)

test/GridapTests/BiharmonicTests.jl

@@ -0,0 +1,54 @@
+module BiharmonicTests
+
+using Test
+using Gridap
+
+# Analytical manufactured solution
+u(x) = x[1]*(x[1]-1)*x[2]*(x[2]-1)
+f(x) = Δ(Δ(u))(x)
+g(x) = Δ(u)(x)
+@test f(VectorValue(0.5,0.5)) == 8.0
+@test g(VectorValue(0.5,0.5)) == -1
+
+# Domain
+domain = (0,1,0,1)
+partition = (4,4)
+model = CartesianDiscreteModel(domain,partition)
+
+# FE space
+order = 2
+V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order),dirichlet_tags="boundary")
+U = TrialFESpace(V,u)
+
+# Triangulation
+Ω = Triangulation(model)
+Γ = BoundaryTriangulation(model)
+Λ = SkeletonTriangulation(model)
+degree = 2*order
+dΩ = Measure(Ω,degree)
+dΓ = Measure(Γ,degree)
+dΛ = Measure(Λ,degree)
+nΓ = get_normal_vector(Γ)
+nΛ = get_normal_vector(Λ)
+
+# Weak form
+const h = (domain[2]-domain[1]) / partition[1]
+const γ = 1
+a(u,v) = ∫( Δ(u)*Δ(v) )dΩ +
+         ∫( - mean(Δ(u))*jump(∇(v)⋅nΛ) - jump(∇(u)⋅nΛ)*mean(Δ(v)) + γ/h*jump(∇(u)⋅nΛ)*jump(∇(v)⋅nΛ) )dΛ
+l(v) = ∫( v*f )dΩ + ∫( g*(∇(v)⋅nΓ) )dΓ
+op = AffineFEOperator(a,l,U,V)
+
+uₕ = solve(op)
+
+# Error
+e = u - uₕ
+l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
+h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))
+el2 = l2(e)
+eh1 = h1(e)
+tol = 1.0e-10
+@test el2 < tol
+@test eh1 < tol
+
+end

test/GridapTests/runtests.jl

@@ -24,4 +24,6 @@ using Test
 
 @time @testset "IsotropicDamage" begin include("IsotropicDamageTests.jl") end
 
+@time @testset "Biharmonic" begin include("BiharmonicTests.jl") end
+
 end # module