Hakai_Project

Package Hakai_Project v0.0.1

This package refactor code from https://github.com/yozoyugen/HAKAI-fem for user to use it.

applyBoundaryConditions!(MODEL, disp_new, step, d_time)

Implements the boundary conditions each time step, using amplitude interpolation.

assignElementMaterial!(PART, INSTANCE, MATERIAL) -> (element_material, element_instance)

Given PART, INSTANCE, and MATERIAL arrays, figure out the material assignments per element and also record which instance each element belongs to.

buildGlobalModel!(PART, INSTANCE) -> (nNode, coordmat, nElement, elementmat)

Given the PART and INSTANCE arrays, iterates over each instance to:

• apply translations/rotations
• assemble global coord and element arrays

Returns the global node count, global coordinate matrix, global element count, and global element matrix.

buildMassDamping(MODEL, elementVolume; mass_scaling=1.0)

Constructs the diagonal mass and damping arrays diagM, diagC based on element densities and volumes. Applies mass scaling.

cal_BVbar_hexa(Pusai_mat, e_position, BVbar)

Given:

• Pusai_mat: a Vector of 8 SMatrix{3,8,Float64} (or MMatrix{3,8}), each storing the local shape function derivatives at one Gauss point of a hexahedron.
• e_position: a 3×8 matrix of the current coordinates of this element’s nodes.
• BVbar: a 6×24 MMatrix (or Array) to which we add the volumetric parts for the B-bar approach.

This function:

1. Loops over each of the 8 Gauss points.
2. Computes the determinant of the Jacobian (detJi) from Pusai_mat[k] * e_position.
3. Accumulates the (1/3)N(detJi) terms into BVbar[1:3, :].
4. Returns the total element volume by summing the 8 detJi values.

The final BVbar is then used in a “Bbar” method (like Simo/Hughes anti-volumetric locking approach).

cal_Bfinal(Bfinal, BVbar, Pusai1, e_position)

Computes the final strain-displacement matrix Bfinal for a single Gauss point of an 8-node hexahedron using a "B-bar" type approach. Specifically:

1. Compute the local 3×3 Jacobian from Pusai1 and e_position.
2. Invert it and build the base B-matrix for the standard part.
3. Subtract or add the volumetric correction from BVbar.
4. Return the determinant of the Jacobian (for weighting).

Arguments: • Bfinal : a 6×24 matrix (initially zero) where we'll store the final B at this Gauss point. • BVbar : a 6×24 matrix containing the volumetric part from cal_BVbar_hexa. • Pusai1 : a 3×8 matrix of shape-function derivatives at this Gauss point (from Pusai_mat). • e_position: a 3×8 matrix of the current coordinates of this element’s 8 nodes.

Returns: • detJi : The determinant of the 3×3 Jacobian, used as the integration weight factor.

Notes: • If you strictly follow your original code logic, you might do: Bfinal = B - BV + BVbar or something similar. Here we embed that logic. • The typical 6×24 B is formed from partial derivatives (dN/dx, etc.) mapped by the Jacobian.

calPusaihexa(integ_num)

Generates the shape-function gradient matrices for an 8-integration-point hexahedral element. Returns an array of size integ_num where each entry is a 3×8 matrix (Pusai1), storing the derivatives of shape functions in local coordinates (gzai, eta, tueta).

integ_num is typically 8 for a standard hexahedron with 2×2×2 integration.

cal_contact_force(c_force3, CT, instance_pair, cp_index,
                  position::Array{Float64,2}, velo::Array{Float64,1}, diag_M::Array{Float64,1},
                  elementMinSize::Float64, elementMaxSize::Float64, d_max::Float64,
                  element_flag::Array{Int,1}, elementmat::Array{Int,2}, bug_report, time_)

Computes contact forces between a set of points (cnodesi) and triangular faces (c_triangles) from the ContactTriangle array CT. The core logic:

1. For each contact pair, gather the relevant nodes/triangles.
2. Build a small coarse grid (nodemapix, nodemapiy, nodemapiz, etc.) to skip far-apart node/triangle checks.
3. For any node that lies within the normal-projection distance d_lim from a triangle, compute normal penalty force and friction. Accumulate into c_force3.
4. We do both damping and friction calculations, so the contact model is basically a penalty approach.

The array c_force3 is 2D: size (fn × Threads.nthreads()) so each thread can accumulate forces without race conditions. The final combination is done after the function, by summation across the second dimension if needed.

cal_node_stress_strain(nNode::Int, elementmat::Array{Int,2}, integ_num::Int, integ_data::IntegDataType)

Averages integration-point (Gauss-point) stress/strain values to each node. Returns a NodeDataType containing:

• node_stress (nNode×6) – the averaged stress at each node • node_strain (nNode×6) – the averaged strain at each node • node_plastic_strain (nNode×6) – (currently unused, but included for completeness) • node_eq_plastic_strain (nNode) – the averaged equivalent plastic strain at each node • node_mises_stress (nNode) – Mises stress computed from node_stress • node_triax_stress (nNode) – triaxial (mean stress / Mises) ratio

Notes:

• integ_data.integ_stress and integ_data.integ_strain are stored in column-major with shape (6 × nElement*integ_num). This function first aggregates them per element, then sums over nodes.

cal_stress_hexa(Qe, integ_stress::Array{Float64,2}, integ_strain::Array{Float64,2},
                integ_yield_stress::Array{Float64,1}, integ_eq_plastic_strain::Array{Float64,1},
                position_::Array{Float64,2}, d_disp_::Array{Float64,1},
                elementmat::Array{Int,2}, element_flag::Array{Int,1}, integ_num::Int,
                Pusai_mat::Vector{SMatrix{3,8,Float64}},  # or your MMatrix type
                MATERIAL::Array{MaterialType,1},
                element_material::Array{Int,1},
                elementMinSize::Float64,
                elementVolume::Array{Float64,1})

Computes the internal forces for each element (Qe) as well as updating integration-point stress/strain with an 8-integration-point hexahedral formulation. The parameters are:

• Qe : a 24×nElement array of element-level force vectors (24 = 3 dofs × 8 nodes). • integ_stress : a 6×(nElementintegnum) array of stress at each Gauss point, updated in-place. • `integstrain` : a 6×(nElementintegnum) array of strain at each Gauss point, updated in-place. • `integyieldstress,integeqplasticstrain: 1D arrays for yield stress and eq. plastic strain, each has nElement*integ_num length. •position: a 3×nNode array of current node coordinates. •ddisp` : a (3×nNode) vector storing this step's incremental displacement (dispnew - dispold). • elementmat: an 8×nElement array listing node indices of each element. • `elementflag: integer array (1 or 0) indicating which elements are active vs. deleted. •integnum: number of integration points per element (8). •Pusaimat: the shape-function derivative data for each integration point. •MATERIAL: array of your customMaterialType, from which we get Dmat, G, plastic table, etc. •element_material: an array giving the material ID for each element. •elementMinSize: minimum element size (used in some checks, can be ignored if not used). •elementVolume` : array of size nElement for (maybe) updated volume if the code modifies it.

What it does:

• Loops over each element (if element_flag[e] == 1).
• Gathers that element’s incremental displacement, forms the B-matrix, and obtains the stress increment.
• If the material is plastic, applies a radial return or linear hardening logic (the code does a simple “incremental” approach).
• Updates integ_stress and integ_strain at each Gauss point.
• Accumulates the element internal force vector into the Qe array.

Returns nothing, but modifies Qe, integ_stress, integ_strain, etc. in-place.

computeElementVolumes(MODEL, Pusai_mat; integ_num=8)

Given the shape function matrix (Pusai_mat) and the MODEL's global coords, compute the volume for each element (summing over integration points). Returns an array of volumes, one per element.

doFractureCheck!(MODEL, element_flag, eq_plastic, triax, integ_stress, integ_strain)

Checks which elements should be "deleted" (elementflag=0) based on the ductile data in MATERIAL. Sets integrstress/strain to zero in any deleted element's integration points. The logic matches your original code.

drawGraph(output)

Plot the output (e.g. a time history) using Plots.

findLineIndices(lines::Vector{String}, pattern::String) -> Vector{Int}

Returns all line indices in lines that contain the given pattern.

get_element_face(MODEL, i)

For the i-th instance in MODEL, returns three arrays describing the faces of each hexahedral element:

• faces : A (nE*6)×4 matrix of node indices, where nE is the number of elements. Each group of 6 rows corresponds to the 6 faces of one element.
• faces_eleid : A vector of length (nE*6) indicating which element each face row belongs to.
• sorted_faces: Same as faces but with each row sorted so that face matching can be performed (useful for detecting shared faces, external faces, etc.).

get_surface_triangle(INSTANCE_i, array_element, contact_element)

Given a single InstanceType (INSTANCE_i), an array of local element indices (array_element), and a subset of those elements (contact_element), returns:

1. c_triangles : An array of triangle faces (two triangles per quadrilateral face).
2. c_triangles_eleid : A matching array of element indices for each triangle row in c_triangles.
3. c_nodes : A sorted and unique list of all node indices used in c_triangles.

This is typically used in contact calculations to extract the exposed surface triangles for a subset of elements, so that point-vs-triangle collision or gap can be computed.

initDisplacementVelocity(MODEL, nNode, d_time, diag_M, mass_scaling)

Creates and returns arrays for disp, dispnew, disppre, velo, and position. Also applies initial conditions from MODEL.IC if present.

main()

Read the input file (from ARGS), run the simulation and produce output.

my3SolveAb(A11::Float64, A21::Float64, A31::Float64,
           A12::Float64, A22::Float64, A32::Float64,
           A13::Float64, A23::Float64, A33::Float64,
           bx::Float64, by::Float64, bz::Float64)

Solves the 3×3 linear system A * x = b, where A is: [ A11 A12 A13 ] [ A21 A22 A23 ] [ A31 A32 A33 ] and b = (bx, by, bz), returning (x1, x2, x3).

It effectively does a manual determinant-based inverse of A:

```math (A^-1) = 1/det(A) * adj(A) then multiplies by b.

my3crossNNz(a1::Float64, a2::Float64, a3::Float64, b1::Float64, b2::Float64, b3::Float64)

Returns (nx, ny, nz) which is the normalized cross product of vectors a = (a1, a2, a3) and b = (b1, b2, b3). In other words, it does:

```math n = (a × b) / ‖a × b‖ where ‖⋅‖ denotes the Euclidean norm. The result (nx, ny, nz) is the unit normal.

Example usage: nx, ny, nz = my3crossNNz(1.0, 0.0, 0.0, 0.0, 1.0, 0.0) # => (0, 0, 1)

parseAmplitudes(lines) -> Vector{AmplitudeType}

Parses *Amplitude blocks.

parseBC(lines, AMPLITUDE, INSTANCE, PART, NSET) -> Vector{BCType}

Parses *Boundary blocks, linking to amplitude if specified.

parseContactFlag(lines) -> Int

Checks if *Contact or *Contact Inclusions (with self-contact) is present. Returns: 0 = no contact 1 = general contact 2 = self-contact

parseContactPairs(lines, SURFACE) -> Vector{CPType}

Parses *Contact Pair blocks into CPType.

parseGlobalElsets(lines, INSTANCE) -> Vector{ELsetType}

Parses any *Elset with "instance=" at the global level.

parseGlobalNsets(lines, INSTANCE) -> Vector{NsetType}

Parses any *Nset with "instance=" at the global level.

parseIC(lines, INSTANCE, PART, NSET) -> Vector{ICType}

Parses *Initial Conditions blocks.

parseInstances(lines, PART) -> Vector{InstanceType}

Finds all *Instance blocks and populates an array of InstanceType.

parseMassScaling(lines) -> Float64

Parses *Fixed Mass Scaling to obtain mass scaling factor.

parseMaterials(lines) -> Vector{MaterialType}

Parses *Material blocks.

parseParts(lines::Vector{String}) -> Vector{PartType}

Parses all *Part sections (with nodes/elements) into an array of PartType.

parseStep(lines) -> (d_time, end_time)

Parses *Dynamic, Explicit parameters.

parseSurfaces(lines, ELSET) -> Vector{SurfaceType}

Parses *Surface blocks and populates SurfaceType.

prepareMaterialProperties!(MODEL)

Loops over MODEL.MATERIAL to compute and store properties like G, Dmat, and detect if there's any fracture. Returns a flag indicating whether fracture is enabled in the model.

readFile(fname::String) -> Vector{String}

Reads the entire input file into a vector of lines.

readInpFile(fname::String) -> ModelType

Top-level function that reads all lines of the .inp file, then parses out all relevant sections into a ModelType.

runMainTimeLoop!(MODEL, diag_M, diag_C, disp, disp_new, disp_pre, velo, position,
                 contactData, elementVolume, Pusai_mat)

Implements the central time-stepping loop:

• updates external forces
• computes contact forces
• enforces boundary conditions
• updates stress/strain
• fracture removal
• writes output periodically

setupContact!(MODEL)

Handles logic to build up contact pairs, surfaces, etc. Returns a struct or NamedTuple with contact data (CT, instancepair, cpindex, etc.).

write_vtk(index::Int, coordmat::AbstractMatrix{T}, elementmat::AbstractMatrix{Int},
          element_flag::AbstractVector{Int}, disp::AbstractVector{T}, velo::AbstractVector{T},
          node_data::NodeDataType)

Writes a .vtk file (in unstructured-grid format) for the mesh and nodal data at the given index. Stores the file in temp/file###.vtk.

Arguments: • coordmat : a 3×nNode matrix of the original node coordinates (column-major). • elementmat : an 8×nElement matrix of node indices for each hexahedral element. • element_flag : an nElement vector of 1/0 flags indicating whether each element is active or deleted. • disp : a 3nNode vector of nodal displacements, stored in [x1,y1,z1, x2,y2,z2, …]. • velo : a 3nNode vector of nodal velocities, same layout as disp. • node_data : a NodeDataType struct, containing node-level stresses, strains, etc.

Output: • Creates a file named "temp/file###.vtk" where ### is the zero-padded index.

Notes: • The VTK file is ASCII, "DATASET UNSTRUCTUREDGRID". • Only hex8 cell types are handled (CELLTYPES 12).