update visualisation and stash Cosserat

tduretz · Dec 7, 2024 · 2bcb3e5 · 2bcb3e5
1 parent 6df84ec
commit 2bcb3e5
Show file tree

Hide file tree

Showing 7 changed files with 1,133 additions and 321 deletions.
diff --git a/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl b/JuliaVisualisation/Main_MatrixCheckSolve_MD7.jl
@@ -42,7 +42,7 @@ function main()
 
     # File
     path     ="/Users/tduretz/REPO/MDOODZ7.0/MDLIB/"
-    filename = string(path, "Stokes_01cpu_step1429_iter00.gzip.h5") 
+    filename = string(path, "Stokes_01cpu_step01_iter00.gzip.h5") 
     # path     ="/Users/tduretz/REPO/MDOODZ7.0/RUNS/RiverTom/"
     # filename = string(path, "Stokes_01cpu_step229_iter00.gzip.h5") 
     # Matrix block A
@@ -75,6 +75,11 @@ function main()
     model = ExtractData(filename,"/model/params");
     nx, nz, Δx, Δz = model
 
+    @show size(MA)
+    @show size(MB)
+    @show size(MC)
+    @show size(MD)
+
     # Spies
     Msc = MA - MB*(MD*MC);
     @show norm(MA-MA')

diff --git a/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ChristmasTreeAnisotropy.jl b/JuliaVisualisation/_SpecificStudies/RiftingAnisotropy/ChristmasTreeAnisotropy.jl
@@ -49,10 +49,10 @@ import Statistics:mean
     # Plasticity
     C          = 5e7
     fric       = 30*π/180
-    dil        = 10*π/180
-    ηvp        = 0.0
+    dil        = 0*π/180
+    ηvp        = 0*1e21
     # Pressure in extension
-    X    = 1 .+ (δ .- 1).*sind.(2*θ) 
+    X    = 1 / (sqrt((δ^2-1)*cos(2*θ)^2 + 1)/δ)
     ϕeff = asin(sin(fric) / X)
     @show ϕeff*180/π
     σH = -Pl*(1-sin(ϕeff))/(1+sin(ϕeff))
@@ -187,15 +187,10 @@ import Statistics:mean
                      2.0*ani*d1 2.0-2.0*ani*d1 2.0*ani*d2;
                     -2.0*ani*d2     2.0*ani*d2 1.0  + 2.0*ani*(d1 - 0.5);]
 
-
-
                 E = [ε̇xxd; ε̇yyd; ε̇xyd]
                 τ0 = [τxx0/ηe; τyy0/ηe; τxy0/ηe]
-                Eeff = E + A\τ0
-
-                @show Eeff
-
-
+                b = A\τ0
+                Eeff = E + b./[1; 1; 2]
 
                 Da11  = 2.0 - 2.0*ani*d1;
                 Da12  = 2.0*ani*d1;
@@ -252,8 +247,8 @@ import Statistics:mean
 
         #####################
 
-        τii_yield     = P*sin(fric) .+ C#*cos(fric)
-        τii_yield_min = τii_yield * (1 / (1 + (δ-1)*sin(2*θ[i]))) 
+        τii_yield     = P*sin(fric) .+ C*cos(fric)
+        τii_yield_min = τii_yield*sqrt((δ^2-1)*cos(2*θ[i])^2 + 1)/δ
 
         ax3 = Axis(f[1, 2], title = L"$$Stress profile ($\delta = 3$)", xlabel = L"$τ_\mathrm{II}$ [GPa]")
         lines!(ax3, τii_rot1./1e9, yc./1e3, label=L"$\sqrt{Y_2^{'}}$", color=:blue )
@@ -268,6 +263,8 @@ import Statistics:mean
         axislegend(position=:rb)
         display(f)
 
+        save("/Users/tduretz/PowerFolders/_manuscripts/RiftingAnisotropy/Figures/ChristmasTreeAnisotropy.png", f, px_per_unit = 4) 
+
         @show τii_rot1./τii_cart1
 
     end