Changeset 26659

Timestamp:

11/23/21 12:28:22 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: getting to Jacobian Determinant

Location:

issm/trunk-jpl/src/jl

Files:

• solve/analyses/stressbalanceanalysis.jl (modified) (1 diff)
• solve/elements.jl (modified) (1 diff)
• solve/gauss.jl (modified) (2 diffs)
• solve/modules.jl (modified) (1 diff)
• solve/nodes.jl (modified) (1 diff)
• solve/solve.jl (modified) (1 diff)
• solve/utils.jl (added)
• solve/vertices.jl (modified) (1 diff)
• test101.jl (modified) (1 diff)

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

@@ -68 +68 @@
 function CreateKMatrix(analysis::StressbalanceAnalysis,element::Tria)# {{{
 
-        #Initialize Element matrix
+        #Internmediaries
+        numnodes = 3
+        
+        #Initialize Element matrix and basis function derivatives
         Ke = ElementMatrix(element.nodes)
+        dbasis = Matrix{Float64}(undef,numnodes,2)
+
+        #Retrieve all inputs and parameters
+        xyz_list = GetVerticesCoordinates(element.vertices)
+        vx_input = GetInput(element, VxEnum)
+        vy_input = GetInput(element, VyEnum)
+        H_input  = GetInput(element, ThicknessEnum)
+
+        #Start integrating
+        gauss = GaussTria(2)
+        for i in 1:gauss.numgauss
+
+                Jdet = JacobianDeterminant(xyz_list, gauss)
+
+                println(Jdet)
+
+                error("STOP")
+                vx = GetInputValue(vx_input, gauss, i)
+                vy = GetInputValue(vy_input, gauss, i)
+        end
+
+
+        error("STOP")
 
         return Ke

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

@@ -102 +102 @@
 
 end # }}}
+
+#Finite Element stuff
+function JacobianDeterminant(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+        #Get Jacobian Matrix
+        J = Jacobian(xyz_list)
+
+        #Get its determinant
+        Jdet = Matrix2x2Determinant(J)
+
+        #check and return
+        if(Jdet<0) error("negative Jacobian Determinant") end
+        return Jdet
+
+end#}}}
+function Jacobian(xyz_list::Matrix{Float64}) #{{{
+
+        J = Matrix{Float64}(undef,2,2)
+
+        x1 = xyz_list[1,1]
+        y1 = xyz_list[1,2]
+        x2 = xyz_list[2,1]
+        y2 = xyz_list[2,2]
+        x3 = xyz_list[3,1]
+        y3 = xyz_list[3,2]
+
+        J[1,1] = .5*(x2-x1)
+        J[1,2] = sqrt(3)/6*(2*x3 -x1 -x2)
+        J[2,1] = .5*(y2-y1)
+        J[2,2] = sqrt(3)/6*(2*y3 -y1 -y2)
+
+        return J
+end#}}}

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

@@ -23 +23 @@
                 coords3 = [0.333333333333333]
         elseif(order==2)
+                npoints = 3
                 weights = [0.577350269189625; 0.577350269189625; 0.577350269189625]
                 coords1 = [0.666666666666667; 0.166666666666667; 0.166666666666667]
@@ -28 +29 @@
                 coords3 = [0.166666666666667; 0.166666666666667; 0.666666666666667]
         elseif(order==3)
-                npoints=4
+                npoints = 4
                 weights = [-0.974278579257493; 0.902109795608790; 0.902109795608790; 0.902109795608790]
                 coords1 = [ 0.333333333333333; 0.600000000000000; 0.200000000000000; 0.200000000000000]

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

@@ -69 +69 @@
         count = 0
         for i in 1:md.mesh.numberofvertices
-                push!(vertices,Vertex(i,x[i],y[i]))
+                push!(vertices,Vertex(i,x[i],y[i],0.))
         end
 

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

@@ -124 +124 @@
                 count = GetDofList(nodes[i],doflist,count,setenum)
         end
-        println(count," ",ndofs)
         @assert count==ndofs
 

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

@@ -14 +14 @@
 include("./modules.jl")
 include("./toolkits.jl")
+include("./utils.jl")
 include("./elementmatrix.jl")
 

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

@@ -4 +4 @@
         x::Float64
         y::Float64
+        z::Float64
 end# }}}
+
+#vertices functions
+function GetVerticesCoordinates(vertices::Vector{Vertex}) #{{{
+
+        #Intermediaries
+        nbv = length(vertices)
+
+        #Allocate
+        xyz_list = Matrix{Float64}(undef,nbv,3)
+
+        #Assign value to xyz_list
+        for i in 1:nbv
+                xyz_list[i,1]=vertices[i].x
+                xyz_list[i,2]=vertices[i].y
+                xyz_list[i,3]=vertices[i].z
+        end
+
+        return xyz_list
+end #}}}

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

@@ -6 +6 @@
 #md = ISSM.triangle(md,"../../test/Exp/Square.exp",180000.)
 md = ISSM.triangle(md,"../../test/Exp/Square.exp",250000.)
+#md = ISSM.triangle(md,"../../test/Exp/Square.exp",1000.)
 md = ISSM.setmask(md,"all","")