dropped number of integration poitns from 12 to 7 in quadratic mortar…

… surfaces intrestingly gives more accurate results, maybe something numerical error in FPG12 integration rule..?
JuliaFEM · Feb 6, 2017 · 40ed4c9 · 40ed4c9
1 parent 9e223c7
commit 40ed4c9
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 17 deletions.
diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl
@@ -623,11 +623,11 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
                     for cell in all_cells
                         virtual_element = Element(Tri3, Int[])
                         update!(virtual_element, "geometry", cell)
-                        for (weight, xi) in get_integration_points(virtual_element, Val{:FPG12})
-                            x_gauss = virtual_element("geometry", xi, time)
+                        for ip in get_integration_points(virtual_element, 3)
+                            x_gauss = virtual_element("geometry", ip, time)
                             xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
-                            detJ = virtual_element(xi, time, Val{:detJ})
-                            w = weight*detJ
+                            detJ = virtual_element(ip, time, Val{:detJ})
+                            w = ip.weight*detJ
                             N1 = vec(slave_element(xi_s, time)*T)
                             De += w*diagm(N1)
                             Me += w*N1*N1'
@@ -728,9 +728,9 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
                     update!(virtual_element, "geometry", cell)
 
                     # 5. loop integration point of integration cell
-                    for (weight, xi) in get_integration_points(virtual_element, Val{:FPG12})
+                    for ip in get_integration_points(virtual_element, 3)
 
-                        x_gauss = virtual_element("geometry", xi, time)
+                        x_gauss = virtual_element("geometry", ip, time)
                         xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
                         xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time)
 
@@ -739,17 +739,19 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
                         N2 = vec(master_element(xi_m, time))
                         Phi = Ae*N1
 
-                        detJ = virtual_element(xi, time, Val{:detJ})
-                        De += weight*Phi*N1'*detJ
-                        Me += weight*Phi*N2'*detJ
+                        detJ = virtual_element(ip, time, Val{:detJ})
+                        w = ip.weight*detJ
+
+                        De += w*Phi*N1'
+                        Me += w*Phi*N2'
                         if props.adjust && haskey(slave_element, "displacement") && haskey(master_element, "displacement")
                             u1 = slave_element("displacement", time)
                             u2 = master_element("displacement", time)
                             xs = N1*(Xs+u1)
                             xm = N2*(Xm+u2)
-                            ge += weight*vec((xm-xs)*Phi')*detJ
+                            ge += w*vec((xm-xs)*Phi')
                         end
-                        area += (weight*detJ)
+                        area += w
                     end # integration points done
 
                 end # integration cells done

diff --git a/test/test_problems_mortar_3d.jl b/test/test_problems_mortar_3d.jl
@@ -144,8 +144,8 @@ end
     maxT = maximum(T)
     stdT = std(T)
     info("minT = $minT, maxT = $maxT, stdT = $stdT")
-    @test maxT - minT < 5.0e-4
-    @test isapprox(stdT, 0.0; atol=1.0e-4)
+    @test maxT - minT < 1.0e-10
+    @test isapprox(stdT, 0.0; atol=1.0e-10)
 
 end
 
@@ -209,7 +209,7 @@ end
     maxabsu3 = maximum(abs(u3))
     stdabsu3 = std(abs(u3))
     info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
-    @test isapprox(stdabsu3, 0.0; atol=1.0e-12)
+    @test isapprox(stdabsu3, 0.0; atol=1.0e-10)
 end
 
 @testset "patch test displacement + abaqus inp + tet10 + adjust + dual basis + alpha=0.2" begin
@@ -274,7 +274,7 @@ end
     maxabsu3 = maximum(abs(u3))
     stdabsu3 = std(abs(u3))
     info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
-    @test isapprox(stdabsu3, 0.0; atol=3.0e-5)
+    @test isapprox(stdabsu3, 0.0; atol=1.0e-6)
 end
 
 @testset "patch test temperature + abaqus inp + tet10 + quadratic surface elements + dual basis + alpha=0.2" begin
@@ -325,6 +325,6 @@ end
     maxT = maximum(T)
     stdT = std(T)
     info("minT = $minT, maxT = $maxT, stdT = $stdT")
-    @test maxT - minT < 5.0e-4
-    @test isapprox(stdT, 0.0; atol=1.0e-4)
+    @test maxT - minT < 1.0e-10
+    @test isapprox(stdT, 0.0; atol=1.0e-10)
 end