Changeset 26675

Timestamp:

11/29/21 07:03:55 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with nodal function derivatives

Location:

issm/trunk-jpl/src/jl/solve

Files:

• analyses/stressbalanceanalysis.jl (modified) (2 diffs)
• elements.jl (modified) (1 diff)
• matice.jl (added)
• solve.jl (modified) (1 diff)
• utils.jl (modified) (1 diff)

issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl

@@ -81 +81 @@
         H_input  = GetInput(element, ThicknessEnum)
 
+        #Prepare material object
+        material = Matice()
+        
         #Start integrating
         gauss = GaussTria(2)
@@ -86 +89 @@
 
                 Jdet = JacobianDeterminant(xyz_list, gauss)
+                NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss)
 
-                println(Jdet)
+                H = GetInputValue(H_input, gauss, i)
 
                 error("STOP")

issm/trunk-jpl/src/jl/solve/elements.jl

@@ -129 +129 @@
 
         J[1,1] = .5*(x2-x1)
-        J[1,2] = sqrt(3)/6*(2*x3 -x1 -x2)
-        J[2,1] = .5*(y2-y1)
+        J[1,2] = .5*(y2-y1)
+        J[2,1] = sqrt(3)/6*(2*x3 -x1 -x2)
         J[2,2] = sqrt(3)/6*(2*y3 -y1 -y2)
 
         return J
 end#}}}
+function JacobianInvert(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+        #Get Jacobian matrix
+        J = Jacobian(xyz_list)
+
+        #Get its determinant
+        Jinv = Matrix2x2Invert(J)
+
+        return Jinv
+end#}}}
+function NodalFunctionsDerivatives(element::Tria,dbasis::Matrix{Float64},xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+        #Get nodal function derivatives in reference element
+        dbasis_ref = Matrix{Float64}(undef,3,2)
+        NodalFunctionsDerivativesReference(dbasis_ref,gauss,P1Enum)
+
+        #Get invert of the Jacobian
+        Jinv = JacobianInvert(xyz_list,gauss)
+
+        #Build dbasis:
+        #[ dNi/dx ] = Jinv * [dNhat_i/dr]
+        #[ dNi/dy ] =        [dNhat_i/ds]
+        for i in 1:3
+                dbasis[i,1] = Jinv[1,1]*dbasis_ref[i,1]+Jinv[1,2]*dbasis_ref[i,2]
+                dbasis[i,2] = Jinv[2,1]*dbasis_ref[i,1]+Jinv[2,2]*dbasis_ref[i,2]
+        end
+
+end#}}}
+function NodalFunctionsDerivativesReference(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{
+
+        if(finiteelement==P0Enum)
+                #Nodal function 1
+                dbasis[1,1]= 0.
+                dbasis[1,2]= 0.
+
+        elseif(finiteelement==P1Enum)
+                #Nodal function 1
+                dbasis[1,1]= -.5
+                dbasis[1,2]= -sqrt(3)/6
+                #Nodal function 2
+                dbasis[2,1]= .5
+                dbasis[2,2]= -sqrt(3)/6
+                #Nodal function 3
+                dbasis[3,1]= 0
+                dbasis[3,2]= sqrt(3)/3
+        else
+                error("Element type ",finiteelement," not supported yet")
+        end
+end#}}}

issm/trunk-jpl/src/jl/solve/solve.jl

@@ -8 +8 @@
 include("./elements.jl")
 include("./constraints.jl")
+include("./matice.jl")
 include("./femmodel.jl")
 include("./analyses/analysis.jl")

issm/trunk-jpl/src/jl/solve/utils.jl

@@ -4 +4 @@
 
 end#}}}
+function Matrix2x2Invert(A::Matrix{Float64}) #{{{
+
+        #Initialize output
+        Ainv = Matrix{Float64}(undef,2,2)
+
+        #Compute determinant
+        det = Matrix2x2Determinant(A)
+        if(abs(det)<eps(Float64)) error("Determinant smaller than machine epsilon") end
+
+        #Multiplication is faster than divsion, so we multiply by the reciprocal
+        det_reciprocal = 1/det
+
+        #compute invert matrix
+        Ainv[1,1]=   A[2,2]*det_reciprocal # =  d/det
+   Ainv[1,2]= - A[1,2]*det_reciprocal # = -b/det
+   Ainv[2,1]= - A[2,1]*det_reciprocal # = -c/det
+   Ainv[2,2]=   A[1,1]*det_reciprocal # =  a/det
+
+        return  Ainv
+
+end#}}}