Problem to create an indexed array - M(1 to 1, 1 to 4){i}

[largurajr]

Hello all, this is my first post here. Well, I need to create a matrix inside a loop. The matrix should have the same name, but different index, for example:
The matrix below, Xc, should be Xc1, Xc2,Xc3,.....XcNE. Now I am storing only the last one, but I need storage all of them.
I don't know how to create an "indexed matrix". Somebody could help me?
While NE <= 200

ReDim Xc(1 To 1, 1 To 4)
Xc(1, 1) = 1 / 8 * N_data(El_data(NE, 2), 2)
Xc(1, 2) = 1 / 8 * N_data(El_data(NE, 3), 2)
Xc(1, 3) = 1 / 8 * N_data(El_data(NE, 4), 2)
Xc(1, 4) = 1 / 8 * N_data(El_data(NE, 5), 2)

Wend

Thank you all in advance!

[SamT]

Redim Xc(1 to 200, 1 to 4)
While NE <= 200
For i = 1 to 4
Xc(NE, i) = N_data(El_data(NE, i + 1), 2) / 8
Next i
Wend

I can't tell what you are trying to do, but what you are doing is

While NE <= 200
For i = 1 to 4
j = El_data(NE, i + 1)
Xc(NE, i) = N_data(j, 2) / 8
Next i
Wend

[snb]

Why do you only provide partial information ?
Is there something to hide ?

[largurajr]

First, thank you for your quickly reply and nope, nothing to hide. It is just because I don't have a code yet, only the beginning. But I am going to send tomorrow.
I am trying to write a finite element code to solve 2D plane stress problems. I don't know if you are familiar with that. So, I need to write the stiffness matrix for each element and then write the global stiffness of the structure. To do that, each element stiffness matrix is K(i), so the total number of matrix depend on the number of elements. I used 200 as an example only. Each element matrix is 8 x 8. Suppose that I have 200 elements then I need K(1), K(2),......K(200). Each one is 8 x 8. Thus, I need write a loop with indexed matrix. What you wrote, is not an indexed matrix. I hope that I clarify better what I have to do. Please, let me know if you understood my question.
Thank you.

[Paul_Hossler]

Nope, didn't understand the question (indexed matrix??) but it sounds like you're looking for an array of arrays(?)

This is a simple 3d version that for demo purposes loads the 200x8x8 elements with random numbers, and then multiplies each the 8x8 to come up with 200 outputs



Option Explicit

Sub Test()
'esentially an array of arrays
Dim K(1 To 200, 1 To 8, 1 To 8) As Double
Dim d1 As Long, d2 As Long, d3 As Long
Dim x As Double

'just to fill in data
For d1 = LBound(K, 1) To UBound(K, 1)
For d2 = LBound(K, 2) To UBound(K, 2)
For d3 = LBound(K, 3) To UBound(K, 3)
K(d1, d2, d3) = 100# * Rnd
Next d3
Next d2
Next d1

'using the data
For d1 = LBound(K, 1) To UBound(K, 1)
x = 0#
For d2 = LBound(K, 2) To UBound(K, 2)
For d3 = LBound(K, 3) To UBound(K, 3)
x = x + K(d1, d2, d3)
Next d3
Next d2

'don't want to do 200
If d1 < 10 Then MsgBox "8x8 product of the " & d1 & "th array is " & x
Next d1
End Sub




And a slightly different way to do it




Sub Test_2()
'esentially an array of arrays
Dim K(1 To 200) As Variant
Dim M(1 To 8, 1 To 8) As Double

Dim d1 As Long, d2 As Long, d3 As Long

'just to fill in data
For d1 = LBound(K) To UBound(K)
Erase M

For d2 = LBound(M, 1) To UBound(M, 1)
For d3 = LBound(M, 1) To UBound(M, 2)
M(d2, d3) = 100# * Rnd
Next d3
Next d2

K(d1) = M
Next d1
End Sub



Now if you were planning to perform matrix math with Excel's functions (e.g. MMULT, MINVERSE, etc.) there's better ways to set it up

Let us know

[mikerickson]

It sounds like you are looking for an array. (K is the name that you used in Post #4.)
Where each element of K is an 8X8 Matrix.


The coding For that would be

Dim k(1 To 200) As Variant
Dim oneMatrix As Variant
Dim i As Long

For i = 1 To 200
ReDim oneMatrix(1 To 8, 1 To 8)
k(i) = oneMatrix
Next i

' put value in
k(2)(1, 4) = 18.32

' get value out
MsgBox k(2)(1, 4)

A more straight forward approach would be to just start off with a 3 dimensional matrix

Dim myArray(1 To 200, 1 To 8, 1 To 8) As Double

' put value in
myArray(2, 1, 4) = 18.32

'get value out
MsgBox myArray(2, 1, 4)

[snb]

A nice tool would be a dictionary, so you can refer to the item's 'Name' (the key).


Sub M_snb()
redim sn(7,7)

with createobject("scripting.dictionary")
for j=1 to 200
.item("Plan_" & j)=sn
next

msgbox .item("Plan_30")(5,5)
end with
End sub

[largurajr]

Thank you.

[largurajr]

Mike, thank you! I think this is what I need.
Thank you all!!!

[mikerickson]

If you need to access the matrices as matrices (e.g. to find their determinant), either the array of arrays technique or snb's dictionary (or collection) would be best.

If you only need to deal with the elements of the matrices, the 3-D array would be easiest.

[largurajr]

Guys, how can I share what I have until now? It is almost nothing, but I would like to share.

[SamT]

Copy the code in the VBA Editor, then Click the # Icon on the VBAExpress Post Editor, then press Ctrl+V.

[largurajr]

Guys, it is just the beginning. I am working on the project. I know the algorithm and I am trying to "translate" to VBA. If you have ideas, please, let me know.


Option Explicit
Sub FEA()
Dim N_data As Variant
Dim El_data As Variant
Dim J1() As Double
Dim J2() As Double
Dim J3() As Double
Dim J4() As Double
Dim Xcoord() As Variant
Dim Xc() As Variant
Dim Ycoord() As Variant
Dim Yc() As Variant
Dim DetJ1() As Variant
Dim DetJ2() As Variant
Dim DetJ3() As Variant
Dim DetJ4() As Variant
Dim XcxJ1() As Variant
Dim XcxJ2() As Variant
Dim XcxJ3() As Variant
Dim XcxJ4() As Variant
Dim BB1() As Variant
Dim BB2() As Variant
Dim BB3() As Variant
Dim BB4() As Variant
Dim BBT1() As Variant
Dim BBT2() As Variant
Dim BBT3() As Variant
Dim BBT4() As Variant
Dim DBJ1() As Variant
Dim BTDB1() As Variant
Dim BTDB2() As Variant
Dim BTDB3() As Variant
Dim BTDB4() As Variant
Dim DB1() As Variant
Dim DB2() As Variant
Dim DB3() As Variant
Dim DB4() As Variant
Dim D() As Variant
Dim Ke() As Variant

Dim A As Integer
Dim b As Integer
Dim i As Integer
Dim ii As Integer
Dim J As Integer
Dim K As Integer
Dim kk As Integer
Dim E As Double
Dim Poisson As Double
Dim NNodes As Double
Dim NElements As Double
Dim NE As Integer
Dim s1 As Double
Dim t1 As Double
Dim s2 As Double
Dim t2 As Double
Dim s3 As Double
Dim t3 As Double
Dim s4 As Double
Dim t4 As Double
Dim a1 As Double
Dim b1 As Double
Dim c1 As Double
Dim d1 As Double
Dim a2 As Double
Dim b2 As Double
Dim c2 As Double
Dim d2 As Double
Dim a3 As Double
Dim b3 As Double
Dim c3 As Double
Dim d3 As Double
Dim a4 As Double
Dim b4 As Double
Dim c4 As Double
Dim d4 As Double
Dim h As Double
K = 5
kk = 4
i = 5
E = Sheets("Input").Cells(1, 3)
Poisson = Sheets("Input").Cells(2, 3)
NNodes = Sheets("Input").Cells(3, 1)
NElements = Sheets("Input").Cells(3, 6)
ReDim Xcoord(NElements)
ReDim Ycoord(NElements)
ReDim XcxJ1(NElements)
ReDim XcxJ2(NElements)
ReDim XcxJ3(NElements)
ReDim XcxJ4(NElements)
ReDim DetJ1(NElements)
ReDim DetJ2(NElements)
ReDim DetJ3(NElements)
ReDim DetJ4(NElements)
ReDim Ke(NElements)
NE = 1
s1 = -0.577350269
t1 = -0.577350269
s2 = -0.577350269
t2 = 0.577350269
s3 = 0.577350269
t3 = 0.577350269
s4 = 0.577350269
t4 = -0.577350269
N_data = Sheets("Input").Range("A5:C" & NNodes + 4).Value 'Nodes data matrix
El_data = Sheets("Input").Range("F5:K" & NElements + 4).Value 'Elements data matrix
'******************************************* MATERIAL BEHAVIOR********************************
If Sheets("Input").Cells(1, 5) = "Plane Stress" Then 'Plane Stress
h = Poisson = Sheets("Input").Cells(1, 8)
ReDim D(1 To 3, 1 To 3)
D(1, 1) = E / (1 - Poisson * Poisson)
D(2, 1) = (Poisson * E) / (1 - Poisson * Poisson)
D(3, 1) = 0
D(1, 2) = (Poisson * E) / (1 - Poisson * Poisson)
D(2, 2) = E / (1 - Poisson * Poisson)
D(3, 2) = 0
D(1, 3) = 0
D(2, 3) = 0
D(3, 3) = ((1 - Poisson) / 2 * E) / (1 - Poisson * Poisson)

Else

h = 1
ReDim D(1 To 3, 1 To 3) 'Plane Strain
D(1, 1) = E / ((1 + Poisson) * (1 - 2 * Poisson)) * (1 - Poisson)
D(2, 1) = Poisson * (E / ((1 + Poisson) * (1 - 2 * Poisson)))
D(3, 1) = 0
D(1, 2) = Poisson * (E / ((1 + Poisson) * (1 - 2 * Poisson)))
D(2, 2) = E / ((1 + Poisson) * (1 - 2 * Poisson)) * (1 - Poisson)
D(3, 2) = 0
D(1, 3) = 0
D(2, 3) = 0
D(3, 3) = E / ((1 + Poisson) * (1 - 2 * Poisson)) * ((1 - 2 * Poisson) / 2)


End If
'******************************************* Element Stiffness matrix *************************************
ReDim J1(1 To 4, 1 To 4)
J1(1, 1) = 0
J1(1, 2) = 1 - t1
J1(1, 3) = t1 - s1
J1(1, 4) = s1 - 1
J1(2, 1) = t1 - 1
J1(2, 2) = 0
J1(2, 3) = s1 + 1
J1(2, 4) = -s1 - t1
J1(3, 1) = s1 - t1
J1(3, 2) = -s1 - 1
J1(3, 3) = 0
J1(3, 4) = t1 + 1
J1(4, 1) = 1 - s1
J1(4, 2) = s1 + t1
J1(4, 3) = -t1 - 1
J1(4, 4) = 0

ReDim J2(1 To 4, 1 To 4)
J2(1, 1) = 0
J2(1, 2) = 1 - t2
J2(1, 3) = t2 - s2
J2(1, 4) = s2 - 1
J2(2, 1) = t2 - 1
J2(2, 2) = 0
J2(2, 3) = s2 + 1
J2(2, 4) = -s2 - t2
J2(3, 1) = s2 - t2
J2(3, 2) = -s2 - 1
J2(3, 3) = 0
J2(3, 4) = t2 + 1
J2(4, 1) = 1 - s2
J2(4, 2) = s2 + t2
J2(4, 3) = -t2 - 1
J2(4, 4) = 0

ReDim J3(1 To 4, 1 To 4)
J3(1, 1) = 0
J3(1, 2) = 1 - t3
J3(1, 3) = t3 - s3
J3(1, 4) = s3 - 1
J3(2, 1) = t3 - 1
J3(2, 2) = 0
J3(2, 3) = s3 + 1
J3(2, 4) = -s3 - t3
J3(3, 1) = s3 - t3
J3(3, 2) = -s3 - 1
J3(3, 3) = 0
J3(3, 4) = t3 + 1
J3(4, 1) = 1 - s3
J3(4, 2) = s3 + t3
J3(4, 3) = -t3 - 1
J3(4, 4) = 0

ReDim J4(1 To 4, 1 To 4)
J4(1, 1) = 0
J4(1, 2) = 1 - t4
J4(1, 3) = t4 - s4
J4(1, 4) = s4 - 1
J4(2, 1) = t4 - 1
J4(2, 2) = 0
J4(2, 3) = s4 + 1
J4(2, 4) = -s4 - t4
J4(3, 1) = s4 - t4
J4(3, 2) = -s4 - 1
J4(3, 3) = 0
J4(3, 4) = t4 + 1
J4(4, 1) = 1 - s4
J4(4, 2) = s4 + t4
J4(4, 3) = -t4 - 1
J4(4, 4) = 0

For NE = 1 To NElements
ReDim Xc(1 To 1, 1 To 4)
Xc(1, 1) = 1 / 8 * N_data(El_data(NE, 2), 2)
Xc(1, 2) = 1 / 8 * N_data(El_data(NE, 3), 2)
Xc(1, 3) = 1 / 8 * N_data(El_data(NE, 4), 2)
Xc(1, 4) = 1 / 8 * N_data(El_data(NE, 5), 2)

Xcoord(NE) = Xc
ReDim Yc(1 To 4, 1 To 1)
Yc(1, 1) = N_data(El_data(NE, 2), 3)
Yc(2, 1) = N_data(El_data(NE, 3), 3)
Yc(3, 1) = N_data(El_data(NE, 4), 3)
Yc(4, 1) = N_data(El_data(NE, 5), 3)

Ycoord(NE) = Yc

XcxJ1 = Application.WorksheetFunction.MMult(Xcoord(NE), J1)

DetJ1 = Application.WorksheetFunction.MMult(XcxJ1, Yc)
XcxJ2 = Application.WorksheetFunction.MMult(Xcoord(NE), J2)

DetJ2 = Application.WorksheetFunction.MMult(XcxJ2, Yc)
XcxJ3 = Application.WorksheetFunction.MMult(Xcoord(NE), J3)

DetJ3 = Application.WorksheetFunction.MMult(XcxJ3, Yc)
XcxJ4 = Application.WorksheetFunction.MMult(Xcoord(NE), J4)

DetJ4 = Application.WorksheetFunction.MMult(XcxJ4, Yc)


a1 = 0.25 * (Yc(1, 1) * (s1 - 1) + Yc(2, 1) * (-1 - s1) + Yc(3, 1) * (1 + s1) + Yc(4, 1) * (1 - s1))
b1 = 0.25 * (Yc(1, 1) * (t1 - 1) + Yc(2, 1) * (1 - t1) + Yc(3, 1) * (1 + t1) + Yc(4, 1) * (-1 - t1))
c1 = 0.25 * (Xc(1, 1) * (t1 - 1) + Xc(1, 2) * (1 - t1) + Xc(1, 3) * (1 + t1) + Xc(1, 4) * (-1 - t1))
d1 = 0.25 * (Xc(1, 1) * (s1 - 1) + Xc(1, 2) * (-1 - s1) + Xc(1, 3) * (1 + s1) + Xc(1, 4) * (1 - s1))

a2 = 0.25 * (Yc(1, 1) * (s2 - 1) + Yc(2, 1) * (-1 - s2) + Yc(3, 1) * (1 + s2) + Yc(4, 1) * (1 - s2))
b2 = 0.25 * (Yc(1, 1) * (t2 - 1) + Yc(2, 1) * (1 - t2) + Yc(3, 1) * (1 + t2) + Yc(4, 1) * (-1 - t2))
c2 = 0.25 * (Xc(1, 1) * (t2 - 1) + Xc(1, 2) * (1 - t2) + Xc(1, 3) * (1 + t2) + Xc(1, 4) * (-1 - t2))
d2 = 0.25 * (Xc(1, 1) * (s2 - 1) + Xc(1, 2) * (-1 - s2) + Xc(1, 3) * (1 + s2) + Xc(1, 4) * (1 - s2))

a3 = 0.25 * (Yc(1, 1) * (s3 - 1) + Yc(2, 1) * (-1 - s3) + Yc(3, 1) * (1 + s3) + Yc(4, 1) * (1 - s3))
b3 = 0.25 * (Yc(1, 1) * (t3 - 1) + Yc(2, 1) * (1 - t3) + Yc(3, 1) * (1 + t3) + Yc(4, 1) * (-1 - t3))
c3 = 0.25 * (Xc(1, 1) * (t3 - 1) + Xc(1, 2) * (1 - t3) + Xc(1, 3) * (1 + t3) + Xc(1, 4) * (-1 - t3))
d3 = 0.25 * (Xc(1, 1) * (s3 - 1) + Xc(1, 2) * (-1 - s3) + Xc(1, 3) * (1 + s3) + Xc(1, 4) * (1 - s3))

a4 = 0.25 * (Yc(1, 1) * (s4 - 1) + Yc(2, 1) * (-1 - s4) + Yc(3, 1) * (1 + s4) + Yc(4, 1) * (1 - s4))
b4 = 0.25 * (Yc(1, 1) * (t4 - 1) + Yc(2, 1) * (1 - t4) + Yc(3, 1) * (1 + t4) + Yc(4, 1) * (-1 - t4))
c4 = 0.25 * (Xc(1, 1) * (t4 - 1) + Xc(1, 2) * (1 - t4) + Xc(1, 3) * (1 + t4) + Xc(1, 4) * (-1 - t4))
d4 = 0.25 * (Xc(1, 1) * (s4 - 1) + Xc(1, 2) * (-1 - s4) + Xc(1, 3) * (1 + s4) + Xc(1, 4) * (1 - s4))

ReDim BB1(1 To 3, 1 To 2)
BB1(1, 1) = (a1 * 0.25 * (t1 - 1) - b1 * 0.25 * (s1 - 1))
BB1(2, 1) = 0
BB1(3, 1) = c1 * 0.25 * (s1 - 1) - d1 * 0.25 * (t1 - 1)
BB1(1, 2) = 0
BB1(2, 2) = c1 * 0.25 * (t1 - 1) - d1 * 0.25 * (s1 - 1)
BB1(3, 2) = a1 * 0.25 * (s1 - 1) - b1 * 0.25 * (t1 - 1)
ReDim BB2(1 To 3, 1 To 2)
BB2(1, 1) = a2 * 0.25 * (t2 - 1) - b2 * 0.25 * (s2 - 1)
BB2(2, 1) = 0
BB2(3, 1) = c2 * 0.25 * (s2 - 1) - d2 * 0.25 * (t2 - 1)
BB2(1, 2) = 0
BB2(2, 2) = c2 * 0.25 * (t2 - 1) - d2 * 0.25 * (s2 - 1)
BB2(3, 2) = a2 * 0.25 * (s2 - 1) - b2 * 0.25 * (t2 - 1)
ReDim BB3(1 To 3, 1 To 2)
BB3(1, 1) = a3 * 0.25 * (t3 - 1) - b3 * 0.25 * (s3 - 1)
BB3(2, 1) = 0
BB3(3, 1) = c3 * 0.25 * (s3 - 1) - d3 * 0.25 * (t3 - 1)
BB3(1, 2) = 0
BB3(2, 2) = c3 * 0.25 * (t3 - 1) - d3 * 0.25 * (s3 - 1)
BB3(3, 2) = a3 * 0.25 * (s3 - 1) - b3 * 0.25 * (t3 - 1)
ReDim BB4(1 To 3, 1 To 2)
BB4(1, 1) = a4 * 0.25 * (t4 - 1) - b4 * 0.25 * (s4 - 1)
BB4(2, 1) = 0
BB4(3, 1) = c4 * 0.25 * (s4 - 1) - d4 * 0.25 * (t4 - 1)
BB4(1, 2) = 0
BB4(2, 2) = c4 * 0.25 * (t4 - 1) - d4 * 0.25 * (s4 - 1)
BB4(3, 2) = a4 * 0.25 * (s4 - 1) - b4 * 0.25 * (t4 - 1)

BBT1 = Application.WorksheetFunction.Transpose(BB1)
BBT2 = Application.WorksheetFunction.Transpose(BB2)
BBT3 = Application.WorksheetFunction.Transpose(BB3)
BBT4 = Application.WorksheetFunction.Transpose(BB4)

DB1 = Application.WorksheetFunction.MMult(D, BB1)
DB2 = Application.WorksheetFunction.MMult(D, BB2)
DB3 = Application.WorksheetFunction.MMult(D, BB3)
DB4 = Application.WorksheetFunction.MMult(D, BB4)

BTDB1 = Application.WorksheetFunction.MMult(BBT1, DB1)
BTDB2 = Application.WorksheetFunction.MMult(BBT2, DB2)
BTDB3 = Application.WorksheetFunction.MMult(BBT3, DB3)
BTDB4 = Application.WorksheetFunction.MMult(BBT4, DB4)

Ke(NE) = Application.WorksheetFunction.Sum(BTDB1, BTDB2, BTDB3, BTDB4)

Next NE

End Sub

[largurajr]

Each K(NE) is the element stiffness. So, I need record the K(1), K(2),......K(NE) because later I am going to write the global stiffness matrix. Each K(NE) is 8x8. I think my code is wrong, I can't retrieve the K(NE) matrix.

[Paul_Hossler]

I know the algorithm and


Can you post it?

It'd probably be easier to follow than code that doesn't seem to be a correct implementation of it

[largurajr]

Basically, the algorithm is below:

15430

The Determinant of the Jacobian matrix is:

[largurajr]

See below.

[largurajr]

The Determinant of the Jacobian Matrix is:
15431

Matrix B is:
15432

Where:

15433

15434
15435

and matrix D is:
15436

Finally Ke is:

15437

Wi = 1 and h = constant
si and ti are constants also.

I need to write the K matrix for each element.

[largurajr]

Guys, I updated the code a little bit:


Sub FEA()
Dim N_data As Variant
Dim El_data As Variant
Dim J1() As Double
Dim J2() As Double
Dim J3() As Double
Dim J4() As Double
Dim Xcoord() As Variant
Dim Xc() As Variant
Dim Ycoord() As Variant
Dim Yc() As Variant
Dim DetJ1 As Double
Dim DetJ2 As Double
Dim DetJ3 As Double
Dim DetJ4 As Double
Dim XcxJ1() As Variant
Dim XcxJ2() As Variant
Dim XcxJ3() As Variant
Dim XcxJ4() As Variant
Dim XcxJ1xYc() As Variant
Dim XcxJ2xYc() As Variant
Dim XcxJ3xYc() As Variant
Dim XcxJ4xYc() As Variant
Dim BB1() As Variant
Dim BB2() As Variant
Dim BB3() As Variant
Dim BB4() As Variant
Dim BBT1() As Variant
Dim BBT2() As Variant
Dim BBT3() As Variant
Dim BBT4() As Variant
Dim DBJ1() As Variant
Dim BTDB1() As Variant
Dim BTDB2() As Variant
Dim BTDB3() As Variant
Dim BTDB4() As Variant
Dim DB1() As Variant
Dim DB2() As Variant
Dim DB3() As Variant
Dim DB4() As Variant
Dim D() As Variant
Dim Ke() As Variant
Dim BTDBJ1() As Variant

Dim A As Integer
Dim b As Integer
Dim i As Integer
Dim ii As Integer
Dim j As Integer
Dim K As Integer
Dim kk As Integer
Dim E As Double
Dim Poisson As Double
Dim NNodes As Double
Dim NElements As Double
Dim NE As Integer
Dim s1 As Double
Dim t1 As Double
Dim s2 As Double
Dim t2 As Double
Dim s3 As Double
Dim t3 As Double
Dim s4 As Double
Dim t4 As Double
Dim a1 As Double
Dim b1 As Double
Dim c1 As Double
Dim d1 As Double
Dim a2 As Double
Dim b2 As Double
Dim c2 As Double
Dim d2 As Double
Dim a3 As Double
Dim b3 As Double
Dim c3 As Double
Dim d3 As Double
Dim a4 As Double
Dim b4 As Double
Dim c4 As Double
Dim d4 As Double
Dim h As Double
K = 5
kk = 4
i = 5
E = Sheets("Input").Cells(1, 3)
Poisson = Sheets("Input").Cells(2, 3)
NNodes = Sheets("Input").Cells(3, 1)
NElements = Sheets("Input").Cells(3, 6)
ReDim Ke(NElements) As Variant
NE = 1
s1 = -0.577350269
t1 = -0.577350269
s2 = -0.577350269
t2 = 0.577350269
s3 = 0.577350269
t3 = -0.577350269
s4 = 0.577350269
t4 = 0.577350269
N_data = Sheets("Input").Range("A5:C" & NNodes + 4).Value 'Nodes data matrix
El_data = Sheets("Input").Range("F5:K" & NElements + 4).Value 'Elements data matrix
'******************************************* MATERIAL BEHAVIOR********************************
If Sheets("Input").Cells(1, 5) = "Plane Stress" Then 'Plane Stress
h = Poisson = Sheets("Input").Cells(1, 8)
ReDim D(1 To 3, 1 To 3)
D(1, 1) = E / (1 - Poisson * Poisson)
D(2, 1) = (Poisson * E) / (1 - Poisson * Poisson)
D(3, 1) = 0
D(1, 2) = (Poisson * E) / (1 - Poisson * Poisson)
D(2, 2) = E / (1 - Poisson * Poisson)
D(3, 2) = 0
D(1, 3) = 0
D(2, 3) = 0
D(3, 3) = ((1 - Poisson) / 2 * E) / (1 - Poisson * Poisson)

Else

h = 1
ReDim D(1 To 3, 1 To 3) 'Plane Strain
D(1, 1) = E / ((1 + Poisson) * (1 - 2 * Poisson)) * (1 - Poisson)
D(2, 1) = Poisson * (E / ((1 + Poisson) * (1 - 2 * Poisson)))
D(3, 1) = 0
D(1, 2) = Poisson * (E / ((1 + Poisson) * (1 - 2 * Poisson)))
D(2, 2) = E / ((1 + Poisson) * (1 - 2 * Poisson)) * (1 - Poisson)
D(3, 2) = 0
D(1, 3) = 0
D(2, 3) = 0
D(3, 3) = E / ((1 + Poisson) * (1 - 2 * Poisson)) * ((1 - 2 * Poisson) / 2)


End If
'******************************************* Element Stiffness matrix *************************************
ReDim J1(1 To 4, 1 To 4)
J1(1, 1) = 0
J1(1, 2) = 1 - t1
J1(1, 3) = t1 - s1
J1(1, 4) = s1 - 1
J1(2, 1) = t1 - 1
J1(2, 2) = 0
J1(2, 3) = s1 + 1
J1(2, 4) = -s1 - t1
J1(3, 1) = s1 - t1
J1(3, 2) = -s1 - 1
J1(3, 3) = 0
J1(3, 4) = t1 + 1
J1(4, 1) = 1 - s1
J1(4, 2) = s1 + t1
J1(4, 3) = -t1 - 1
J1(4, 4) = 0

ReDim J2(1 To 4, 1 To 4)
J2(1, 1) = 0
J2(1, 2) = 1 - t2
J2(1, 3) = t2 - s2
J2(1, 4) = s2 - 1
J2(2, 1) = t2 - 1
J2(2, 2) = 0
J2(2, 3) = s2 + 1
J2(2, 4) = -s2 - t2
J2(3, 1) = s2 - t2
J2(3, 2) = -s2 - 1
J2(3, 3) = 0
J2(3, 4) = t2 + 1
J2(4, 1) = 1 - s2
J2(4, 2) = s2 + t2
J2(4, 3) = -t2 - 1
J2(4, 4) = 0

ReDim J3(1 To 4, 1 To 4)
J3(1, 1) = 0
J3(1, 2) = 1 - t3
J3(1, 3) = t3 - s3
J3(1, 4) = s3 - 1
J3(2, 1) = t3 - 1
J3(2, 2) = 0
J3(2, 3) = s3 + 1
J3(2, 4) = -s3 - t3
J3(3, 1) = s3 - t3
J3(3, 2) = -s3 - 1
J3(3, 3) = 0
J3(3, 4) = t3 + 1
J3(4, 1) = 1 - s3
J3(4, 2) = s3 + t3
J3(4, 3) = -t3 - 1
J3(4, 4) = 0

ReDim J4(1 To 4, 1 To 4)
J4(1, 1) = 0
J4(1, 2) = 1 - t4
J4(1, 3) = t4 - s4
J4(1, 4) = s4 - 1
J4(2, 1) = t4 - 1
J4(2, 2) = 0
J4(2, 3) = s4 + 1
J4(2, 4) = -s4 - t4
J4(3, 1) = s4 - t4
J4(3, 2) = -s4 - 1
J4(3, 3) = 0
J4(3, 4) = t4 + 1
J4(4, 1) = 1 - s4
J4(4, 2) = s4 + t4
J4(4, 3) = -t4 - 1
J4(4, 4) = 0

For NE = 1 To NElements
ReDim Xc(1 To 1, 1 To 4)
Xc(1, 1) = N_data(El_data(NE, 2), 2)
Xc(1, 2) = N_data(El_data(NE, 3), 2)
Xc(1, 3) = N_data(El_data(NE, 4), 2)
Xc(1, 4) = N_data(El_data(NE, 5), 2)

ReDim Xcdet(1 To 1, 1 To 4) 'coordinates divided by 8 to calculate the determinant of J
Xcdet(1, 1) = 1 / 8 * N_data(El_data(NE, 2), 2)
Xcdet(1, 2) = 1 / 8 * N_data(El_data(NE, 3), 2)
Xcdet(1, 3) = 1 / 8 * N_data(El_data(NE, 4), 2)
Xcdet(1, 4) = 1 / 8 * N_data(El_data(NE, 5), 2)

ReDim Yc(1 To 4, 1 To 1)
Yc(1, 1) = N_data(El_data(NE, 2), 3)
Yc(2, 1) = N_data(El_data(NE, 3), 3)
Yc(3, 1) = N_data(El_data(NE, 4), 3)
Yc(4, 1) = N_data(El_data(NE, 5), 3)


XcxJ1 = Application.WorksheetFunction.MMult(Xcdet, J1)

XcxJ1xYc = Application.WorksheetFunction.MMult(XcxJ1, Yc)
XcxJ2 = Application.WorksheetFunction.MMult(Xcdet, J2)

XcxJ2xYc = Application.WorksheetFunction.MMult(XcxJ2, Yc)
XcxJ3 = Application.WorksheetFunction.MMult(Xcdet, J3)

XcxJ3xYc = Application.WorksheetFunction.MMult(XcxJ3, Yc)
XcxJ4 = Application.WorksheetFunction.MMult(Xcdet, J4)

XcxJ4xYc = Application.WorksheetFunction.MMult(XcxJ4, Yc)

DetJ1 = XcxJ1xYc(1)
DetJ2 = XcxJ2xYc(1)
DetJ3 = XcxJ3xYc(1)
DetJ4 = XcxJ4xYc(1)


a1 = 0.25 * (Yc(1, 1) * (s1 - 1) + Yc(2, 1) * (-1 - s1) + Yc(3, 1) * (1 + s1) + Yc(4, 1) * (1 - s1))
b1 = 0.25 * (Yc(1, 1) * (t1 - 1) + Yc(2, 1) * (1 - t1) + Yc(3, 1) * (1 + t1) + Yc(4, 1) * (-1 - t1))
c1 = 0.25 * (Xc(1, 1) * (t1 - 1) + Xc(1, 2) * (1 - t1) + Xc(1, 3) * (1 + t1) + Xc(1, 4) * (-1 - t1))
d1 = 0.25 * (Xc(1, 1) * (s1 - 1) + Xc(1, 2) * (-1 - s1) + Xc(1, 3) * (1 + s1) + Xc(1, 4) * (1 - s1))

a2 = 0.25 * (Yc(1, 1) * (s2 - 1) + Yc(2, 1) * (-1 - s2) + Yc(3, 1) * (1 + s2) + Yc(4, 1) * (1 - s2))
b2 = 0.25 * (Yc(1, 1) * (t2 - 1) + Yc(2, 1) * (1 - t2) + Yc(3, 1) * (1 + t2) + Yc(4, 1) * (-1 - t2))
c2 = 0.25 * (Xc(1, 1) * (t2 - 1) + Xc(1, 2) * (1 - t2) + Xc(1, 3) * (1 + t2) + Xc(1, 4) * (-1 - t2))
d2 = 0.25 * (Xc(1, 1) * (s2 - 1) + Xc(1, 2) * (-1 - s2) + Xc(1, 3) * (1 + s2) + Xc(1, 4) * (1 - s2))

a3 = 0.25 * (Yc(1, 1) * (s3 - 1) + Yc(2, 1) * (-1 - s3) + Yc(3, 1) * (1 + s3) + Yc(4, 1) * (1 - s3))
b3 = 0.25 * (Yc(1, 1) * (t3 - 1) + Yc(2, 1) * (1 - t3) + Yc(3, 1) * (1 + t3) + Yc(4, 1) * (-1 - t3))
c3 = 0.25 * (Xc(1, 1) * (t3 - 1) + Xc(1, 2) * (1 - t3) + Xc(1, 3) * (1 + t3) + Xc(1, 4) * (-1 - t3))
d3 = 0.25 * (Xc(1, 1) * (s3 - 1) + Xc(1, 2) * (-1 - s3) + Xc(1, 3) * (1 + s3) + Xc(1, 4) * (1 - s3))

a4 = 0.25 * (Yc(1, 1) * (s4 - 1) + Yc(2, 1) * (-1 - s4) + Yc(3, 1) * (1 + s4) + Yc(4, 1) * (1 - s4))
b4 = 0.25 * (Yc(1, 1) * (t4 - 1) + Yc(2, 1) * (1 - t4) + Yc(3, 1) * (1 + t4) + Yc(4, 1) * (-1 - t4))
c4 = 0.25 * (Xc(1, 1) * (t4 - 1) + Xc(1, 2) * (1 - t4) + Xc(1, 3) * (1 + t4) + Xc(1, 4) * (-1 - t4))
d4 = 0.25 * (Xc(1, 1) * (s4 - 1) + Xc(1, 2) * (-1 - s4) + Xc(1, 3) * (1 + s4) + Xc(1, 4) * (1 - s4))

ReDim BB1(1 To 3, 1 To 8)
BB1(1, 1) = ((a1 * 0.25 * (t1 - 1) - b1 * 0.25 * (s1 - 1))) / DetJ1
BB1(2, 1) = 0
BB1(3, 1) = (c1 * 0.25 * (s1 - 1) - d1 * 0.25 * (t1 - 1)) / DetJ1
BB1(1, 2) = 0
BB1(2, 2) = (c1 * 0.25 * (s1 - 1) - d1 * 0.25 * (s1 - 1)) / DetJ1
BB1(3, 2) = (a1 * 0.25 * (t1 - 1) - b1 * 0.25 * (t1 - 1)) / DetJ1
BB1(1, 3) = (a2 * 0.25 * (1 - t1) - b2 * 0.25 * (-1 - s1)) / DetJ1
BB1(2, 3) = 0
BB1(3, 3) = (c2 * 0.25 * (-1 - s1) - d2 * 0.25 * (1 - t1)) / DetJ1
BB1(1, 4) = 0
BB1(2, 4) = (c2 * 0.25 * (-1 - s1) - d2 * 0.25 * (1 - t1)) / DetJ1
BB1(3, 4) = (a2 * 0.25 * (1 - t1) - b2 * 0.25 * (-1 - s1)) / DetJ1
BB1(1, 5) = (a3 * 0.25 * (1 + t1) - b3 * 0.25 * (1 + s1)) / DetJ1
BB1(2, 5) = 0
BB1(3, 5) = (c3 * 0.25 * (1 + s1) - d3 * 0.25 * (1 + t1)) / DetJ1
BB1(1, 6) = 0
BB1(2, 6) = (c3 * 0.25 * (1 + s1) - d3 * 0.25 * (1 + t1)) / DetJ1
BB1(3, 6) = (a3 * 0.25 * (1 + t1) - b3 * 0.25 * (1 + s1)) / DetJ1
BB1(1, 7) = (a4 * 0.25 * (-1 - t1) - b4 * 0.25 * (1 - s1)) / DetJ1
BB1(2, 7) = 0
BB1(3, 7) = (c4 * 0.25 * (1 - s1) - d4 * 0.25 * (-1 - t1)) / DetJ1
BB1(1, 8) = 0
BB1(2, 8) = (c4 * 0.25 * (1 - s1) - d4 * 0.25 * (-1 - t1)) / DetJ1
BB1(3, 8) = (a4 * 0.25 * (-1 - t1) - b4 * 0.25 * (1 - s1)) / DetJ1


ReDim BB2(1 To 3, 1 To 8)
BB2(1, 1) = ((a1 * 0.25 * (t2 - 1) - b1 * 0.25 * (s2 - 1))) / DetJ2
BB2(2, 1) = 0
BB2(3, 1) = (c1 * 0.25 * (s2 - 1) - d1 * 0.25 * (t2 - 1)) / DetJ2
BB2(1, 2) = 0
BB2(2, 2) = (c1 * 0.25 * (s2 - 1) - d1 * 0.25 * (s2 - 1)) / DetJ2
BB2(3, 2) = (a1 * 0.25 * (t2 - 1) - b1 * 0.25 * (t2 - 1)) / DetJ2
BB2(1, 3) = (a2 * 0.25 * (1 - t2) - b2 * 0.25 * (-1 - s2)) / DetJ2
BB2(2, 3) = 0
BB2(3, 3) = (c2 * 0.25 * (-1 - s2) - d2 * 0.25 * (1 - t2)) / DetJ2
BB2(1, 4) = 0
BB2(2, 4) = (c2 * 0.25 * (-1 - s2) - d2 * 0.25 * (1 - t2)) / DetJ2
BB2(3, 4) = (a2 * 0.25 * (1 - t2) - b2 * 0.25 * (-1 - s2)) / DetJ2
BB2(1, 5) = (a3 * 0.25 * (1 + t2) - b3 * 0.25 * (1 + s2)) / DetJ2
BB2(2, 5) = 0
BB2(3, 5) = (c3 * 0.25 * (1 + s2) - d3 * 0.25 * (1 + t2)) / DetJ2
BB2(1, 6) = 0
BB2(2, 6) = (c3 * 0.25 * (1 + s2) - d3 * 0.25 * (1 + t2)) / DetJ2
BB2(3, 6) = (a3 * 0.25 * (1 + t2) - b3 * 0.25 * (1 + s2)) / DetJ2
BB2(1, 7) = (a4 * 0.25 * (-1 - t2) - b4 * 0.25 * (1 - s2)) / DetJ2
BB2(2, 7) = 0
BB2(3, 7) = (c4 * 0.25 * (1 - s2) - d4 * 0.25 * (-1 - t2)) / DetJ2
BB2(1, 8) = 0
BB2(2, 8) = (c4 * 0.25 * (1 - s2) - d4 * 0.25 * (-1 - t2)) / DetJ2
BB2(3, 8) = (a4 * 0.25 * (-1 - t2) - b4 * 0.25 * (1 - s2)) / DetJ2


ReDim BB3(1 To 3, 1 To 8)
BB3(1, 1) = ((a1 * 0.25 * (t3 - 1) - b1 * 0.25 * (s3 - 1))) / DetJ3
BB3(2, 1) = 0
BB3(3, 1) = (c1 * 0.25 * (s3 - 1) - d1 * 0.25 * (t3 - 1)) / DetJ3
BB3(1, 2) = 0
BB3(2, 2) = (c1 * 0.25 * (s3 - 1) - d1 * 0.25 * (s3 - 1)) / DetJ3
BB3(3, 2) = (a1 * 0.25 * (t3 - 1) - b1 * 0.25 * (t3 - 1)) / DetJ3
BB3(1, 3) = (a2 * 0.25 * (1 - t3) - b2 * 0.25 * (-1 - s3)) / DetJ3
BB3(2, 3) = 0
BB3(3, 3) = (c2 * 0.25 * (-1 - s3) - d2 * 0.25 * (1 - t3)) / DetJ3
BB3(1, 4) = 0
BB3(2, 4) = (c2 * 0.25 * (-1 - s3) - d2 * 0.25 * (1 - t3)) / DetJ3
BB3(3, 4) = (a2 * 0.25 * (1 - t3) - b2 * 0.25 * (-1 - s3)) / DetJ3
BB3(1, 5) = (a3 * 0.25 * (1 + t3) - b3 * 0.25 * (1 + s3)) / DetJ3
BB3(2, 5) = 0
BB3(3, 5) = (c3 * 0.25 * (1 + s3) - d3 * 0.25 * (1 + t3)) / DetJ3
BB3(1, 6) = 0
BB3(2, 6) = (c3 * 0.25 * (1 + s3) - d3 * 0.25 * (1 + t3)) / DetJ3
BB3(3, 6) = (a3 * 0.25 * (1 + t3) - b3 * 0.25 * (1 + s3)) / DetJ3
BB3(1, 7) = (a4 * 0.25 * (-1 - t3) - b4 * 0.25 * (1 - s3)) / DetJ3
BB3(2, 7) = 0
BB3(3, 7) = (c4 * 0.25 * (1 - s3) - d4 * 0.25 * (-1 - t3)) / DetJ3
BB3(1, 8) = 0
BB3(2, 8) = (c4 * 0.25 * (1 - s3) - d4 * 0.25 * (-1 - t3)) / DetJ3
BB3(3, 8) = (a4 * 0.25 * (-1 - t3) - b4 * 0.25 * (1 - s3)) / DetJ3


ReDim BB4(1 To 3, 1 To 8)
BB4(1, 1) = ((a1 * 0.25 * (t4 - 1) - b1 * 0.25 * (s4 - 1))) / DetJ4
BB4(2, 1) = 0
BB4(3, 1) = (c1 * 0.25 * (s4 - 1) - d1 * 0.25 * (t4 - 1)) / DetJ4
BB4(1, 2) = 0
BB4(2, 2) = (c1 * 0.25 * (s4 - 1) - d1 * 0.25 * (s4 - 1)) / DetJ4
BB4(3, 2) = (a1 * 0.25 * (t4 - 1) - b1 * 0.25 * (t4 - 1)) / DetJ4
BB4(1, 3) = (a2 * 0.25 * (1 - t4) - b2 * 0.25 * (-1 - s4)) / DetJ4
BB4(2, 3) = 0
BB4(3, 3) = (c2 * 0.25 * (-1 - s4) - d2 * 0.25 * (1 - t4)) / DetJ4
BB4(1, 4) = 0
BB4(2, 4) = (c2 * 0.25 * (-1 - s4) - d2 * 0.25 * (1 - t4)) / DetJ4
BB4(3, 4) = (a2 * 0.25 * (1 - t4) - b2 * 0.25 * (-1 - s4)) / DetJ4
BB4(1, 5) = (a3 * 0.25 * (1 + t4) - b3 * 0.25 * (1 + s4)) / DetJ4
BB4(2, 5) = 0
BB4(3, 5) = (c3 * 0.25 * (1 + s4) - d3 * 0.25 * (1 + t4)) / DetJ4
BB4(1, 6) = 0
BB4(2, 6) = (c3 * 0.25 * (1 + s4) - d3 * 0.25 * (1 + t4)) / DetJ4
BB4(3, 6) = (a3 * 0.25 * (1 + t4) - b3 * 0.25 * (1 + s4)) / DetJ4
BB4(1, 7) = (a4 * 0.25 * (-1 - t4) - b4 * 0.25 * (1 - s4)) / DetJ4
BB4(2, 7) = 0
BB4(3, 7) = (c4 * 0.25 * (1 - s4) - d4 * 0.25 * (-1 - t4)) / DetJ4
BB4(1, 8) = 0
BB4(2, 8) = (c4 * 0.25 * (1 - s4) - d4 * 0.25 * (-1 - t4)) / DetJ4
BB4(3, 8) = (a4 * 0.25 * (-1 - t4) - b4 * 0.25 * (1 - s4)) / DetJ4

BBT1 = Application.WorksheetFunction.Transpose(BB1)
BBT2 = Application.WorksheetFunction.Transpose(BB2)
BBT3 = Application.WorksheetFunction.Transpose(BB3)
BBT4 = Application.WorksheetFunction.Transpose(BB4)

DB1 = Application.WorksheetFunction.MMult(D, BB1)
DB2 = Application.WorksheetFunction.MMult(D, BB2)
DB3 = Application.WorksheetFunction.MMult(D, BB3)
DB4 = Application.WorksheetFunction.MMult(D, BB4)

BTDB1 = Application.WorksheetFunction.MMult(BBT1, DB1)
BTDB2 = Application.WorksheetFunction.MMult(BBT2, DB2)
BTDB3 = Application.WorksheetFunction.MMult(BBT3, DB3)
BTDB4 = Application.WorksheetFunction.MMult(BBT4, DB4)

Ke(NE) = Application.WorksheetFunction.Sum(BTDB1, BTDB2, BTDB3, BTDB4)

Next NE

End Sub

[Paul_Hossler]

A couple of quick observations

1. Ke appears to be intended as 4 element matrix, but

Ke(NE) = Application.WorksheetFunction.Sum(BTDB1, BTDB2, BTDB3, BTDB4)

is not how to do matrix addition. Each element needs to be added

(MMULT can can make a 3x3 matrix in a variant)


You can step through the code Sub Right below if you want


2. Instead of

Dim BB1() As Variant
ReDim BB1(1 To 3, 1 To 2)


I think it'd be clearer to use

Dim BB1(1 to 3, 1 to 2) As Double


3. You have a lot of 'four-ples' in your structure




Dim J1() As Double
Dim J2() As Double
Dim J3() As Double
Dim J4() As Double
Dim DetJ1() As Variant
Dim DetJ2() As Variant
Dim DetJ3() As Variant
Dim DetJ4() As Variant
Dim XcxJ1() As Variant
Dim XcxJ2() As Variant
Dim XcxJ3() As Variant
Dim XcxJ4() As Variant
Dim BB1() As Variant
Dim BB2() As Variant
Dim BB3() As Variant
Dim BB4() As Variant
Dim BBT1() As Variant
Dim BBT2() As Variant
Dim BBT3() As Variant
Dim BBT4() As Variant
Dim DBJ1() As Variant
Dim BTDB1() As Variant
Dim BTDB2() As Variant
Dim BTDB3() As Variant
Dim BTDB4() As Variant
Dim DB1() As Variant
Dim DB2() As Variant
Dim DB3() As Variant
Dim DB4() As Variant



you could simplify your logic a little more with

Dim DB(1 to 4, 1 to 3, 1 to 3) As Double or depending
Dim DB(1 to 4, 1 to 3) As Double

if you wanted to




Option Explicit
Sub Wrong()
Dim NE As Long
Dim Ke() As Variant
Dim BTDB1(1 To 3) As Long
Dim BTDB2(1 To 3) As Long
Dim BTDB3(1 To 3) As Long
Dim BTDB4(1 To 3) As Long

Dim i As Long, j As Long

NE = 3
ReDim Ke(1 To NE)

For j = 1 To NE

'fill data
For i = 1 To 3
BTDB1(i) = 11
BTDB2(i) = 11
BTDB3(i) = 11
BTDB4(i) = 11
Next I

'Ke is a Variant array (1..3) with a scaler (double = 132)
Ke(j) = Application.WorksheetFunction.Sum(BTDB1, BTDB2, BTDB3, BTDB4)

Next j

Stop
End Sub

Sub Right()
Dim NE As Long
Dim Ke() As Variant
Dim BTDB1(1 To 3) As Long
Dim BTDB2(1 To 3) As Long
Dim BTDB3(1 To 3) As Long
Dim BTDB4(1 To 3) As Long

Dim i As Long, j As Long

NE = 3
ReDim Ke(1 To NE)

For j = 1 To NE

'fill data
For i = 1 To 3
BTDB1(i) = 11
BTDB2(i) = 11
BTDB3(i) = 11
BTDB4(i) = 11
Next I

'Ke is a Variant array with 3 elements, each element being a 1 dim array, 1..3
Ke(j) = BTDB1
'now that each Ke is an array, we can add the others
For i = 1 To 3
Ke(j)(i) = Ke(j)(i) + BTDB2(i)
Ke(j)(i) = Ke(j)(i) + BTDB3(i)
Ke(j)(i) = Ke(j)(i) + BTDB4(i)
Next i
Next j

Stop

End Sub


Sub Right2()
Dim BTDB1(1 To 3, 1 To 3) As Long
Dim BTDB2(1 To 3, 1 To 3) As Long

Dim V As Variant

Dim i As Long, j As Long

'fill data
For i = 1 To 3
For j = 1 To 3
BTDB1(i, j) = 100 * i + j
BTDB2(i, j) = 100 * i + j
Next j
Next I

'variant with 3x3 matrix
V = Application.WorksheetFunction.MMult(BTDB1, BTDB2)
Stop

End Sub

[SamT]

I think tat a much easier way is to use Excel as much as possible. To that end I have started a workbook demonstrating what I mean.

When this workbook is opened, Calculation in Excel will be set to Manual only. Click the "Calculate" button on sheet "Math" to make Excel Calculate, When the workbook is closed, Excels Calculation is restored to Automatic.

I have Defined Names so that the Formulas closely resemble the formulas in your code. Note that some Names are suffixed with an underscore. This is required to create a valid Name. These Names can be used in Code. Ex: Dim J1 As Variant: Set J1 = Range("J1_"), creates an array, J1, with the values in Range("J1_")

15438