Merge pull request #604 from gridap/improve_divergence

Implemented divergence for rank 2 tensor-valued functions

Author: fverdugo

src/Fields/AutoDiff.jl

@@ -68,7 +68,49 @@ function symmetric_gradient(f::Function,x::Point)
 end
 
 function divergence(f::Function,x::Point)
-  tr(gradient(f,x))
+  divergence(f,x,return_value(f,x))
+end
+
+function divergence(f::Function,x::Point,fx)
+  tr(gradient(f,x,fx))
+end
+
+function divergence(f::Function,x::Point,fx::TensorValue{2,2})
+  g(x) = SVector(f(x).data)
+  a = ForwardDiff.jacobian(g,get_array(x))
+  VectorValue(
+    a[1,1]+a[2,2],
+    a[3,1]+a[4,2],
+  )
+end
+
+function divergence(f::Function,x::Point,fx::TensorValue{3,3})
+  g(x) = SVector(f(x).data)
+  a = ForwardDiff.jacobian(g,get_array(x))
+  VectorValue(
+    a[1,1]+a[2,2]+a[3,3],
+    a[4,1]+a[5,2]+a[6,3],
+    a[7,1]+a[8,2]+a[9,3],
+   )
+end
+
+function divergence(f::Function,x::Point,fx::SymTensorValue{2})
+  g(x) = SVector(f(x).data)
+  a = ForwardDiff.jacobian(g,get_array(x))
+  VectorValue(
+    a[1,1]+a[2,2],
+    a[2,1]+a[3,2],
+  )
+end
+
+function divergence(f::Function,x::Point,fx::SymTensorValue{3})
+  g(x) = SVector(f(x).data)
+  a = ForwardDiff.jacobian(g,get_array(x))
+  VectorValue(
+    a[1,1]+a[2,2]+a[3,3],
+    a[2,1]+a[4,2]+a[5,3],
+    a[3,1]+a[5,2]+a[6,3],
+   )
 end
 
 function curl(f::Function,x::Point)

test/FieldsTests/DiffOperatorsTests.jl

@@ -93,6 +93,9 @@ for x in xs
   @test (∇×u)(x) == grad2curl(∇u(x))
   @test Δ(u)(x) == Δu(x)
   @test ε(u)(x) == εu(x)
+  @test ∇(u)(x) == ∇u(x)
+  @test Δ(u)(x) == (∇⋅∇u)(x)
+  @test (∇⋅εu)(x) == VectorValue(4.,-1.)
 end
 
 u(x) = VectorValue( x[1]^2 + 2*x[2]^2, 0 )
@@ -105,6 +108,24 @@ for x in xs
   @test (∇×u)(x) == grad2curl(∇u(x))
   @test Δ(u)(x) == Δu(x)
   @test ε(u)(x) == εu(x)
+  @test Δ(u)(x) == (∇⋅∇u)(x)
+  @test (∇⋅εu)(x) == VectorValue(4.,0.)
+end
+
+u(x) = VectorValue( x[1]^2 + 2*x[3]^2, -x[1]^2, -x[2]^2 + x[3]^2 )
+∇u(x) = TensorValue( 2*x[1],0,4*x[3],  -2*x[1],0,0,  0,-2*x[2],2*x[3] )
+Δu(x) = VectorValue( 6, -2, 0 )
+εu(x) = symmetric_part(∇u(x))
+
+xs = [ Point(1.,1.,2.0), Point(2.,0.,1.), Point(0.,3.,0.), Point(-1.,3.,2.)]
+for x in xs
+  @test ∇(u)(x) == ∇u(x)
+  @test (∇⋅u)(x) == tr(∇u(x)) 
+  @test (∇×u)(x) == grad2curl(∇u(x))
+  @test Δ(u)(x) == Δu(x)
+  @test ε(u)(x) == εu(x)
+  @test Δ(u)(x) == (∇⋅∇u)(x)
+  @test (∇⋅εu)(x) == VectorValue(4.,-1.,1.)
 end
 
 end # module