Fixed volume calculation for Abaqus

CHANGELOG.md

@@ -8,10 +8,16 @@ SPDX-License-Identifier: BSD-3-Clause
 
 All notable changes to this project will be documented in this file.
 
-## [1.1.3] - 2024-
+## [1.1.3] - 2024-03-28
 
 ### Fixed
 
+- Volume calculation for Abaqus
+- #137
+
+### Added 
+
+- Abaqus Test
 
 ## [1.1.2] - 2024-03-25
 

Project.toml

@@ -9,6 +9,7 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+Delaunay = "07eb4e4e-0c6d-46ef-bc4e-83d5e5d860a9"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Exodus = "f57ae99e-f805-4780-bdca-96e224be1e5a"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
@@ -41,6 +42,7 @@ ArgParse = "^1"
 CSV = "^0.10"
 DataFrames = "^1"
 Dates = "^1"
+Delaunay = "^1"
 Dierckx = "^0.5"
 Exodus = "^0.12"
 FastGaussQuadrature = "^1"

src/IO/mesh_data.jl

@@ -5,6 +5,7 @@
 using LinearAlgebra
 using AbaqusReader
 using DataFrames
+using Delaunay
 using NearestNeighbors: BallTree
 using NearestNeighbors: inrange
 using PrettyTables
@@ -526,10 +527,17 @@ function read_mesh(filename::String, params::Dict)
         nodes = mesh["nodes"]
         elements = mesh["elements"]
         element_sets = mesh["element_sets"]
+        element_types = mesh["element_types"]
 
-        dof = size(nodes[1])[1]
         dof = 2
 
+        nodes_vector = collect(values(nodes))
+        node_z = [node[end] for node in nodes_vector]
+        if length(unique(node_z)) > 2
+            dof = 3
+        end
+        @info "Abaqus mesh with $dof DOF"
+
         mesh_df = ifelse(dof == 2,
             DataFrame(x=[], y=[], volume=[], block_id=Int[]),
             DataFrame(x=[], y=[], z=[], volume=[], block_id=Int[])
@@ -541,7 +549,8 @@ function read_mesh(filename::String, params::Dict)
         nsets = Dict{String,Any}()
 
         # sort element_sets by length
-        element_sets_keys = sort(collect(keys(element_sets)), by=x -> length(element_sets[x]), rev=true)
+        # element_sets_keys = sort(collect(keys(element_sets)), by=x -> length(element_sets[x]), rev=true)
+        element_sets_keys = collect(keys(element_sets))
         for key in element_sets_keys
             ns_nodes = Int64[]
             for (i, element_id) in enumerate(element_sets[key])
@@ -551,13 +560,13 @@ function read_mesh(filename::String, params::Dict)
                 end
                 push!(ns_nodes, id)
                 node_ids = elements[element_id]
+                element_type = element_types[element_id]
                 vertices = [nodes[node_id] for node_id in node_ids]
+                volume = calculate_volume(string(element_type), vertices)
                 center = sum(vertices) / size(vertices)[1]
                 if dof == 2
-                    volume = area_of_polygon(vertices)
                     push!(mesh_df, (x=center[1], y=center[2], volume=volume, block_id=block_id))
                 else
-                    volume = abs(dot(vertices[1] - vertices[4], cross(vertices[2] - vertices[4], vertices[3] - vertices[4]))) / 6.0
                     push!(mesh_df, (x=center[1], y=center[2], z=center[3], volume=volume, block_id=block_id))
                 end
                 push!(element_written, element_id)
@@ -586,6 +595,50 @@ function read_mesh(filename::String, params::Dict)
 
 end
 
+"""
+    calculate_volume(element_type::String, vertices::Vector{Vector{Float64}})
+
+Calculate the volume of a element.
+
+# Arguments
+- `element_type`: The element type of the element.
+- `vertices`: The vertices of the element.
+# Returns
+- `volume`: The volume of the element.
+"""
+function calculate_volume(element_type::String, vertices::Vector{Vector{Float64}})
+    if element_type == "Quad4"
+        return area_of_polygon(vertices)
+    elseif element_type == "Tet4"
+        return tetrahedron_volume(vertices)
+    elseif element_type == "Hex8"
+        # Convert vertices to a 3x8 matrix for easier manipulation
+        vertices_matrix = hcat(vertices...)
+        tetrahedra_indices = delaunay(Matrix(vertices_matrix'))
+        volumes = [tetrahedron_volume(tetrahedra_indices.points[tetrahedra_indices.simplices[i, :], :]) for i in 1:size(tetrahedra_indices.simplices)[1]]
+        return sum(volumes)
+    else
+        @error "Element type $element_type currently not supported"
+        return nothing
+    end
+end
+
+"""
+    tetrahedron_volume(vertices)
+
+Calculate the volume of a tetrahedron.
+
+# Arguments
+- `vertices`: The vertices of the tetrahedron.
+# Returns
+- `volume`: The volume of the tetrahedron.
+"""
+function tetrahedron_volume(vertices)
+    mat = hcat(vertices, ones(4))  # Augmenting matrix with ones in the fourth column
+    volume = abs(det(mat) / 6)   # Using det function to calculate determinant
+    return volume
+end
+
 """
     area_of_polygon(vertices)
 

test/fullscale_tests/test_Abaqus/Abaqus.yaml

@@ -4,27 +4,27 @@
 
 PeriLab:
   Discretization:
-    Input Mesh File: Dogbone.inp
+    Input Mesh File: plate.inp
     Type: Abaqus
     Distribution Type: Neighbor based
   Physics:
     Material Models:
       Material_1:
         Material Model: PD Solid Elastic
-        Young's Modulus: 200000.0
+        Young's Modulus: 100000.0
         Poisson's Ratio: 0.34
   Blocks:
     block_1:
       Material Model: Material_1
-      Horizon: 8
+      Horizon: 4
       Density: 2.7e-9
     block_2:
       Material Model: Material_1
-      Horizon: 8
+      Horizon: 4
       Density: 2.7e-9
     block_3:
       Material Model: Material_1
-      Horizon: 8
+      Horizon: 4
       Density: 2.7e-9
   Boundary Conditions:
     Displacement-1:
@@ -41,7 +41,7 @@ PeriLab:
     Output1:
       Output File Type: Exodus
       Output Filename: Abaqus
-      Number of Output Steps: 100
+      Number of Output Steps: 200
       Output Variables:
         Displacements: True
         Number of Neighbors: True