-
Notifications
You must be signed in to change notification settings - Fork 0
/
ProblemDefinition.m
66 lines (54 loc) · 1.43 KB
/
ProblemDefinition.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
function ProblemDefinition()
global nDim nNodes nElements nNodesElement nDoF nEquations ...
nEdgesElement Coord ID IEN LM EBC Params f;
%Solution paramaters
nDoF=3;
nDim=3;
nEdgesElement=6; %this is actually the number of faces in 3D
%Mesh parameters
Params.Nx = 2^10; % Number of elements along x-axis
Params.Ny = 2^1; % Number of elements along y-axis
Params.Nz = 2^1; % Number of elements along z-axis
% Parameters for Plate Geometry
Params.L = 2^7; % Length of plate, i.e. (0 < x < L)
Params.c = 1; % Half-height of plate, i.e. (-c < y < c)
Params.t = 1; % Thickness of plate (0<z<t)
% Parameters for Material Properties
Params.E = 10^7; % Youngs Modulus
Params.v = 0.25; % Poisson Ratio
% Parameters for Loading Conditions
Params.P = 80; % Shear load at beam end
%Mesh paramaters
Nx = Params.Nx;
Ny = Params.Ny;
Nz = Params.Nz;
L = Params.L;
c = Params.c;
t= Params.t;
% Build Mesh
x = linspace(0,L,Nx+1);
y = linspace(-c,c,Ny+1);
z = linspace(-t,t,Nz+1);
[X,Y,Z] = meshgrid(x,y,z);
X = permute(X,[2 1 3]);
Y = permute(Y,[2 1 3]);
Z=permute(Z,[2 1 3]);
Coord = [X(:), Y(:), Z(:)];
% Build IEN Array
IEN=linIEN(Nx,Ny,Nz);
nNodes=size(Coord,1);
nElements=size(IEN,2);
nNodesElement=size(IEN,1);
init_data();
%Allocate arrays
ID=zeros(nDoF,nNodes);
LM=zeros(nNodesElement*nDoF,nElements);
%Create ID
I=full(EBC'==0);
nEquations=sum(sum(I));
ID(I)=1:nEquations;
ID=ID';
%Create LM
P=ID(IEN,:)';
LM(:)=P(:);
end