Materials and Computational Mechanics

What is XFEM?

2013-05-30T01:30:00.001-05:00

XFEM, Extended Finite Element Method, is a method used to model discontinuity within a continuum with the aid of additional degrees of freedom, it is specifically used to model crack propagation throughout the model. It connects LEFM to finite element by adding 8 additional DOFs to the nodes of the element containing the crack tip, and the jump in displacement across the cracked surface is modeled by adding heaviside function to the nodes of the elements containing the crack.

There are a few free XFEM codes such as XFEM++, getfem++ and openxfem++.

It is also implemented in Abaqus.

It is also possible to connect EPFM with FE by using J integral, but this is only valid for local plasticity.

Isotropic Hardening, von Mises plasticity, UMAT

2013-02-12T19:39:00.001-06:00

UMAT by

Terminal

2012-05-16T21:31:00.000-05:00

After Installing Ubuntu, you are mainly going to work through Terminal, therefore I will discuss about basic commands in Terminal.
First of all, in order to open a new terminal, you can go to dash home at the top of the side bar at the left of your screen and type terminal and then click on it or you can click Ctrl+Alt+T
After you opened a terminal, you can check the current working directory by typing pwd.
In order to see what files and folders are in this directory type ls. Remember that there might be some options for the command that you are typing, you can read about these options through manual pages, e.g. type man ls. You are able to move up and down in the manual by Up and Down arrows of your keyboard and you can search by typing the favorable text after slash e.g. /full. You can quit the manual by typing q.

installing Ubuntu

2012-05-16T20:44:00.002-05:00

In order to install Ubuntu, please visit this website
After you are directed to the website, click on Ubuntu Desktop, and it will get you to the download page for the latest version of Ubuntu (12.04) which I believe is still the Beta version.
I believe that it is better to download Ubuntu 10.4 instead of 12.04.
Anyways, just follow the instructions on the Ubuntu web page, and remember that it is possible to have Ubuntu alongside your other OS.

ubuntu

2012-05-16T20:31:00.000-05:00

I am going to write a few posts about Ubuntu, which is a free open source operating system and there are a wide variety of open source applications and softwares available for this OS, e.g. GNU compilers, tcl tk. With the help of these softwares we can provide ourselves with more comprehensive tools for programming.
I start with beginner lessons and move forward to more advanced topics towards programming in C, fortran, tcl tk, Octave and etc.

Structural Reliability course Summary Part 3

2012-03-30T17:39:00.001-05:00

SR-summary-part3 (click to download from Scribd)

Structural Reliability course Summary Part 2

2012-03-30T14:44:00.000-05:00

SR-Summary-Part2 (click to download from Scribd)

Structural Reliability course summary part1

2012-03-30T14:37:00.001-05:00

Hello all,

I am uploading the structural reliability course summary part by part, I hope you find it useful.
SR-Summary-Part1 (click to download the file from Scribd)

FEniCS

2012-01-16T22:18:00.003-06:00

Hi again and Happy New Year to you all,

In this post I am introducing the FEniCS software which is indeed an open source software.

As a brief introduction: "The FEniCS Project is a collaborative project for the development of innovative concepts and tools for automated scientific computing, with a particular focus on automated solution of differential equations by finite element methods."

This project is made possible by different contributors worldwide. Contributors

You can use notepad ++ to code and there are two options to interact with FEniCS: python and c++.

There are a few manuals and books on the web, from which I am using Tutorials

You can download the software from http://fenicsproject.org/

programming

2011-12-09T17:49:00.001-06:00

Hi Guys!

I found an interesting app on google chrome webstore which contains lots and lots of step by step lessons about programming languages, specifically C and C++. Its fascinating to know that there are videos, pdf files, cpp files and examples available on the website too. So if you are interested you can visit the website at:
website
The website describes itself as "WiBit.net is a video tutorial web site offering cutting edge programming and computer tutorials."
It is good to know that they use free compiler GCC and notepad++ as the text editor.

J2-Hardening plasticity

2011-12-03T02:37:00.001-06:00

Below you may find the constitutive driver for the J2-Plasticity model with isotropic hardening, considering associated flow rule, which means that the flow potential is the same as the yield surface and therefore the normal to the yield surface is identical to the normal to the plastic potential surface.

function [sig,eps,state_new,E_tang]=J2hardening(oldsig,oldeps,Deps,state_old,modpar)
% function [sig,eps,state_new,E_tang]=J2hardening(oldsig,oldeps,Deps,state_old,modpar)
%----------- Description ---------
% J2 elasto plasticity with isotropic hardening constitutive driver
% ------------ INPUT -------------
% oldsig = stress vector from the last increment [9X1]
% oldeps = strain vector from the last increment [9X1]
% Deps = strain increment [9X1]
% state_old = structure array with the specified fields and values
% could contain eps_e, eps_p, kappa, mu, etc.
% modpar = contains material constants
% ----------- OUTPUT -------------
% sig = stress vector from the current increment [9X1]
% eps = strain vector from the current increment [9X1]
% state_new = updated structured values
% E_tang = Tangent stiffness matrix
% Written by
% Shah, 02 December 2011

E=modpar.E;
nu=modpar.nu;
H=modpar.H;
sy=modpar.sy;
oldeps_p=state_old.eps_p;
oldkap=state_old.kap;
K=E/3/(1-2*nu);
G=E/2/(1+nu);
v1=[1 1 1 0 0 0 0 0 0]';
v2=[0 0 0 1 1 1 1 1 1]';
I=diag(v1+0.5*v2);
Idev=I-(1/3)*(v1*v1');
Ee=2*G*Idev+K*(v1*v1');
I2=diag(v1+v2);
I2dev=I2-(1/3)*(v1*v1');
sigT=oldsig+Ee*Deps;
sigDT=I2dev*sigT;
sigET=sqrt(1.5)*norm(sigDT);
phiT=sigET-sy-oldkap;
sig_m1=(sigT(1)+sigT(2)+sigT(3))/3;
eps=oldeps+Deps;
if phiT<=0
sig=sigT;
mu=0;
kap=oldkap;
eps_p=oldeps_p;
E_tang=Ee;
else
mu=phiT/(3*G+H);
c=1-mu*3*G/sigET;
sigD=sigDT*c;
sig_m=(sigT(1)+sigT(2)+sigT(3))/3;
sig=sigD+sig_m.*v1;
nuT=(1.5)*(sigDT/sigET);
kap=oldkap+H*mu;
E_tang=Ee-2*G*((2*G/(3*G+H)).*((oldkap+sy)/sigET)*nuT*nuT'+mu*3*G*Idev/sigET);
eps_p=oldeps_p+mu;
end
state_new.eps_p=eps_p;
state_new.J2=0;
state_new.I1=0;
state_new.kap=kap;
state_new.mu=mu;

GID

2011-12-02T18:59:00.001-06:00

There is another mesh generator GID which is quite a powerful and complete mesh tool. This software is not free, but you can download it for free and use it with some limitations such as number of elements or nodes. Here I post the download link, give it a try, you won't be disappointed.
http://www.gidhome.com/download
There are versions for Windows, Mac and Linux for 32bit and 64bit.

Transfer your mesh from GMSH to Matlab

2011-11-18T01:55:00.001-06:00

I believe examples are on of the best ways to learn, so I start from an example.

1- build your geometry in GMSH (you can easily do that by studying the first tutorial of the software)

2- Mesh the geometry by going to the Mesh tab and clicking on the 2D. This will build a triangular non structured mesh

3- Extract the mesh by going to File-Save Mesh, or you can use the hot key (Shift+Ctrl+S). This will give you a file with the same name as the project name and with the extension "MSH"

4- Open the file with Wordpad, Notepad or etc.

$MeshFormat

2.2 0 8

$EndMeshFormat

$Nodes

1 0 0 0

2 4 0 0

3 4 4 0

4 0 4 0

5 2 2 0

6 2.5 2 0

7 1.5 2 0

8 1.333333333330429 0 0

9 2.666666666663595 0 0

10 4 1.333333333330429 0

11 4 2.666666666663595 0

12 2.666666666669443 4 0

13 1.33333333333796 4 0

14 0 2.666666666669443 0

15 0 1.33333333333796 0

16 2.43301270189563 2.249999999994092 0

17 2.250000000005625 2.433012701888972 0

18 1.999999999999019 2.5 0

19 1.749999999993242 2.433012701888317 0

20 1.566987298104043 2.249999999993525 0

21 1.566987298103185 1.75000000000796 0

22 1.749999999991474 1.566987298112704 0

23 1.999999999998037 1.5 0

24 2.25000000000626 1.566987298111395 0

25 2.433012701896161 1.750000000006827 0

26 0.9761017850604321 2.346814881605136 0

27 3.02389821493962 1.65318511839525 0

28 2.334500907280026 3.134455102784824 0

29 1.665499092719324 0.8655448972157629 0

30 2.839249308207946 2.495807560436933 0

31 1.160750691791375 1.504192439564072 0

32 1.510100448723277 2.957936003268833 0

33 2.489899551276439 1.04206399673127 0

34 3.168083376431483 3.259385865977071 0

35 0.7639071134243337 3.194283510308682 0

36 3.236092886575931 0.8057164896913898 0

37 0.8319166235682255 0.7406141340235588 0

$EndNodes

$Elements

1 15 2 0 1 1

2 15 2 0 2 2

3 15 2 0 3 3

4 15 2 0 4 4

5 15 2 0 5 5

6 15 2 0 6 6

7 15 2 0 7 7

8 1 2 0 1 1 8

9 1 2 0 1 8 9

10 1 2 0 1 9 2

11 1 2 0 2 2 10

12 1 2 0 2 10 11

13 1 2 0 2 11 3

14 1 2 0 3 3 12

15 1 2 0 3 12 13

16 1 2 0 3 13 4

17 1 2 0 4 4 14

18 1 2 0 4 14 15

19 1 2 0 4 15 1

20 1 2 0 5 6 16

21 1 2 0 5 16 17

22 1 2 0 5 17 18

23 1 2 0 5 18 19

24 1 2 0 5 19 20

25 1 2 0 5 20 7

26 1 2 0 6 7 21

27 1 2 0 6 21 22

28 1 2 0 6 22 23

29 1 2 0 6 23 24

30 1 2 0 6 24 25

31 1 2 0 6 25 6

32 2 2 0 9 11 3 34

33 2 2 0 9 34 3 12

34 2 2 0 9 13 4 35

35 2 2 0 9 35 4 14

36 2 2 0 9 9 2 36

37 2 2 0 9 36 2 10

38 2 2 0 9 37 1 8

39 2 2 0 9 15 1 37

40 2 2 0 9 32 13 35

41 2 2 0 9 33 9 36

42 2 2 0 9 34 12 28

43 2 2 0 9 37 8 29

44 2 2 0 9 31 21 7

45 2 2 0 9 6 25 27

46 2 2 0 9 10 11 27

47 2 2 0 9 30 16 6

48 2 2 0 9 6 27 30

49 2 2 0 9 30 27 11

50 2 2 0 9 7 20 26

51 2 2 0 9 7 26 31

52 2 2 0 9 14 15 26

53 2 2 0 9 31 26 15

54 2 2 0 9 37 31 15

55 2 2 0 9 22 21 31

56 2 2 0 9 22 31 29

57 2 2 0 9 31 37 29

58 2 2 0 9 17 16 30

59 2 2 0 9 17 30 28

60 2 2 0 9 34 30 11

61 2 2 0 9 30 34 28

62 2 2 0 9 14 26 35

63 2 2 0 9 32 35 26

64 2 2 0 9 10 27 36

65 2 2 0 9 33 36 27

66 2 2 0 9 32 19 18

67 2 2 0 9 8 9 29

68 2 2 0 9 23 22 29

69 2 2 0 9 18 17 28

70 2 2 0 9 18 28 32

71 2 2 0 9 12 13 28

72 2 2 0 9 32 28 13

73 2 2 0 9 33 29 9

74 2 2 0 9 33 24 23

75 2 2 0 9 23 29 33

76 2 2 0 9 20 19 32

77 2 2 0 9 20 32 26

78 2 2 0 9 25 24 33

79 2 2 0 9 25 33 27

$EndElements

5- You can see the nodes from the line $Nodes to $EndNodes. The first number below $Nodes is the number of the nodes, The first number at the left column specifies the node number and the rest are the x,y and z position of the nodes (The nodes are highlighted in yellow). Copy the nodes into a separate "txt" file and name is as Nodes.txt.

6- The elements start from $Elements to $EndElements. In this example we are only interested in 2D elements, these elements start from the line in which the second column's value is 2. 32 2 2 0 9 11 3 34. Copy the elements into a separate "txt" file and name it as Elements.txt. The corresponding nodes of the element are the last three numbers, e.g. 32 2 2 0 9 11 3 3 4.

7- Copy these two files into the working directory of your Matlab code.

8- Use the following code to transfer the data into Matlab:

----------------------------------------------------------

%----------- Description ---------

% Transfer Mesh from GMSH to Matlab

%---------------------------------

% Written by

% Shah, 18 November 2011

clc

clear all

close all

load Elements.txt

load Nodes.txt

nrelm=size(Elements,1);

nnod=size(Nodes,1);

Edof=zeros(nrelm,4);

Ex=zeros(nrelm,3); Ey=Ex; Ez=Ex;

Edof(:,1)=sort(1:nrelm);

Elm(:,1)=sort(1:nrelm);

for i=1:nrelm

Edof(i,2:7)=[2*Elements(i,6)-1 Elements(i,6)*2 2*Elements(i,7)-1 ...

Elements(i,7)*2 2*Elements(i,8)-1 2*Elements(i,8)];

Elm(i,2:4)=Elements(i,6:8);

Ex(i,:)=[Nodes(Elm(i,2),2) Nodes(Elm(i,3),2) Nodes(Elm(i,4),2)];

Ey(i,:)=[Nodes(Elm(i,2),3) Nodes(Elm(i,3),3) Nodes(Elm(i,4),3)];

Ez(i,:)=[Nodes(Elm(i,2),4) Nodes(Elm(i,3),4) Nodes(Elm(i,4),4)];

end

nrdof=max(max(Edof));

figure(1);

hold on;

for i=1:length(Edof(:,1));

plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));

end

axis equal

------------------------------------------------------------------------------------------------

The output in Matlab looks like this:

This does not look a like good approximation of the geometry, a finer mesh would look like:

Homework 8 solution

2011-11-17T23:02:00.001-06:00

Homework 8 solution (click to download from Scribd)

HW#8

2011-11-17T22:58:00.001-06:00

HW8 problem description (click to download from Scribd)

GMSH

2011-11-13T15:38:00.001-06:00

Here is the link to a very useful free mesh generator and post processor software (GMSH) which is distributed under the terms of the GNU General Public License (GPL). http://geuz.org/gmsh/ There are versions available for Win, Mac and Linux.
You can find tutorials in the tutorial folder of the downloaded file.

Below is an example of unstructured mesh for a circular surface with a squared hole inside.

4 node quadrilateral element for non symmetrical rectangles

2011-11-11T03:52:00.001-06:00

Below you may find the Matlab function for 4 node quadrilateral element which is capable of meshing non symmetric rectangles. I will try to do the same for triangles.

function [Edof,Ex,Ey]=quadmesh2(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=quadmesh2(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates quadrilateral 4 node elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. elements in x direction
% ny= No. elements in y direction
% 4 3
% o----------o
% | |
% | |
% | |
% o-------- o
% 1 2
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 4]
% Ey= element y coordinates matrix [nrelem X 4]
% Edof= element DOF [nrelem X 9]
% Written by
% Shah, 11 November 2011

ld=abs(xcorner(2)-xcorner(1));
lu=abs(xcorner(3)-xcorner(4));
hr=abs(ycorner(3)-ycorner(2));
hl=abs(ycorner(4)-ycorner(1));
for i=1:ny
for j=1:nx
dl=(i-1)*(nx+1)+j;
ul=dl+(nx+1);
dr=dl+1;
ur=ul+1;
l1=(i-1)*(lu-ld)/ny+ld; l2=i*(lu-ld)/ny+ld;
h1=(j-1)*(hr-hl)/nx+hl; h2=j*(hr-hl)/nx+hl;
x_cd=(i-1)*(xcorner(4)-xcorner(1))/ny+xcorner(1);
x_cu=i*(xcorner(4)-xcorner(1))/ny+xcorner(1);
y_cl=(j-1)*(ycorner(2)-ycorner(1))/nx+ycorner(1);
y_cr=j*(ycorner(2)-ycorner(1))/nx+ycorner(1);
xy_dl=[(j-1)*(l1/nx)+x_cd (i-1)*(h1/ny)+y_cl];
xy_dr=[j*(l1/nx)+x_cd (i-1)*(h2/ny)+y_cr];
xy_ul=[(j-1)*(l2/nx)+x_cu i*(h1/ny)+y_cl];
xy_ur=[j*(l2/nx)+x_cu i*(h2/ny)+y_cr];
el=(i-1)*nx+j;
Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];
Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ur(1) xy_ul(1)];
Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ur(2) xy_ul(2)];
end
end

figure(1);
hold on;
for i=1:length(Edof(:,1));
plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));
end
axis equal

3 node symmetric triangular elements

2011-11-11T01:59:00.001-06:00

Below you may find a Matlab function to generate 3 node symmetric triangular elements as I promised.
Note that nx and ny are the number of rectangular elements in each direction.

function [Edof,Ex,Ey]=symtrimesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=symtrimesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates triangular 3 node symmetric elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. of rectangulat elements in x direction
% ny= No. of rectangulat elements in y direction
% 7 8 9
% o---------o---------o
% | + | + |
% | + | + |
% | + | + |
% | + | + |
% 4o---------5---------o6
% | + | + |
% | + | + |
% | + | + |
% | + | + |
% o---------o---------o
% 1 2 3
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 3]
% Ey= element y coordinates matrix [nrelem X 3]
% Edof= element DOF [nrelem X 7]
% Written by
% Shah, 11 November 2011

l=xcorner(2)-xcorner(1);

h=ycorner(3)-ycorner(2);

for i=1:ny
for j=1:nx

dl=(i-1)*(2*nx+1)+2*j-1;
dm=dl+1;
dr=dl+2;
ul=dl+(2*nx+1);
um=ul+1;
ur=ul+2;
xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_dm=(xy_dl+xy_dr)./2;
xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
xy_um=(xy_ul+xy_ur)./2;
el=4*nx*(i-1)+4*(j-1)+1;
if floor(i/2)-i/2~=0
Edof(el,:)=[el,2*dl-1,2*dl,2*um-1,2*um,2*ul-1,2*ul];
Edof(el+1,:)=[el+1,2*dl-1,2*dl,2*dm-1,2*dm,2*um-1,2*um];
Edof(el+2,:)=[el+2,2*dm-1,2*dm,2*dr-1,2*dr,2*um-1,2*um];
Edof(el+3,:)=[el+3,2*dr-1,2*dr,2*ur-1,2*ur,2*um-1,2*um];
Ex(el,:)=[xy_dl(1) xy_um(1) xy_ul(1)];
Ex(el+1,:)=[xy_dl(1) xy_dm(1) xy_um(1)];
Ex(el+2,:)=[xy_dm(1) xy_dr(1) xy_um(1)];
Ex(el+3,:)=[xy_dr(1) xy_ur(1) xy_um(1)];
Ey(el,:)=[xy_dl(2) xy_um(2) xy_ul(2)];
Ey(el+1,:)=[xy_dl(2) xy_dm(2) xy_um(2)];
Ey(el+2,:)=[xy_dm(2) xy_dr(2) xy_um(2)];
Ey(el+3,:)=[xy_dr(2) xy_ur(2) xy_um(2)];
else
Edof(el,:)=[el,2*dl-1,2*dl,2*dm-1,2*dm,2*ul-1,2*ul];
Edof(el+1,:)=[el+1,2*dm-1,2*dm,2*um-1,2*um,2*ul-1,2*ul];
Edof(el+2,:)=[el+2,2*dm-1,2*dm,2*ur-1,2*ur,2*um-1,2*um];
Edof(el+3,:)=[el+3,2*dm-1,2*dm,2*dr-1,2*dr,2*ur-1,2*ur];
Ex(el,:)=[xy_dl(1) xy_dm(1) xy_ul(1)];
Ex(el+1,:)=[xy_dm(1) xy_um(1) xy_ul(1)];
Ex(el+2,:)=[xy_dm(1) xy_ur(1) xy_um(1)];
Ex(el+3,:)=[xy_dm(1) xy_dr(1) xy_ur(1)];
Ey(el,:)=[xy_dl(2) xy_dm(2) xy_ul(2)];
Ey(el+1,:)=[xy_dm(2) xy_um(2) xy_ul(2)];
Ey(el+2,:)=[xy_dm(2) xy_ur(2) xy_um(2)];
Ey(el+3,:)=[xy_dm(2) xy_dr(2) xy_ur(2)];

end
end
end

figure(1);

hold on;
for i=1:length(Edof(:,1));

plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));

end
axis equal

3 node triangular element

2011-11-11T00:28:00.001-06:00

Below you may find a Matlab function to generate 3 node triangular elements. The output will not be symmetric, so I will try to write one for symmetric triangular elements soon.
Note that nx and ny are the number of rectangular elements in each direction.

function [Edof,Ex,Ey]=trimesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=trimesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates triangular 3 node elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. elements in x direction
% ny= No. elements in y direction
% 3 4
% o---------o
% | + |
% | + |
% | + |
% | + |
% o---------o
% 1 2
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 3]
% Ey= element y coordinates matrix [nrelem X 3]
% Edof= element DOF [nrelem X 7]
% Written by
% Shah, 11 November 2011

l=xcorner(2)-xcorner(1);

h=ycorner(3)-ycorner(2);

for i=1:ny
for j=1:nx

dl=(i-1)*(nx+1)+j;

ul=dl+(nx+1);

dr=dl+1;

ur=ul+1;

xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];

xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];

xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];

xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];

el=2*(i-1)*nx+2*j-1;

Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ul-1,2*ul];

Edof(el+1,:)=[el+1,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];

Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ul(1)];

Ex(el+1,:)=[xy_dr(1) xy_ur(1) xy_ul(1)];

Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ul(2)];

Ey(el+1,:)=[xy_dr(2) xy_ur(2) xy_ul(2)];

end
end

figure(1);

hold on;
for i=1:length(Edof(:,1));

plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));

end
axis equal

4 node quadrilateral element

2011-11-10T01:08:00.003-06:00

Below you may find a Matlab function written by me to generate 4 node quadrilateral elements. It is very easy and straight forward.

function [Edof,Ex,Ey]=quadmesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=quadmesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates quadrilateral 4 node elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left]
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. elements in x direction
% ny= No. elements in y direction
% 4 3
% o---------o
% | |
% | |
% | |
% o---------o
% 1 2
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 4]
% Ey= element y coordinates matrix [nrelem X 4]
% Edof= element DOF [nrelem X 9]
% Written by
% Shah, 10 November 2011

l=xcorner(2)-xcorner(1);
h=ycorner(3)-ycorner(2);
for i=1:ny
for j=1:nx
dl=(i-1)*(nx+1)+j;
ul=dl+(nx+1);
dr=dl+1;
ur=ul+1;
xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
el=(i-1)*nx+j;
Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];
Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ur(1) xy_ul(1)];
Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ur(2) xy_ul(2)];
end
end

figure(1);
hold on;
for i=1:length(Edof(:,1));
plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));
end
axis equal

Homework 7 solution

2011-11-09T23:11:00.000-06:00

Homework 7 solution (click to download from Scribd)

HW#7

2011-11-09T22:58:00.002-06:00

HW7 problem description (click to download from Scribd)

Homework 6 solution

2011-10-28T03:32:00.001-05:00

Homework 6 solution (click to download from Scribd)

HW#6

2011-10-28T03:30:00.003-05:00

HW6 problem description (click to download from Scribd)

Homework 5 solution

2011-10-28T03:28:00.002-05:00

Homework 5 solution (click to download from Scribd)