Example:Example Quad4i 01
From Nemesis
(Difference between revisions)
m |
|||
| (3 intermediate revisions by one user not shown) | |||
| Line 2: | Line 2: | ||
# /****************************************************************************** | # /****************************************************************************** | ||
# * nemesis. an experimental finite element code. * | # * nemesis. an experimental finite element code. * | ||
| − | # * Copyright (C) 2004- | + | # * Copyright (C) 2004-2008 F.E.Karaoulanis [http://www.nemesis-project.org] * |
# * * | # * * | ||
# * This program is free software; you can redistribute it and/or modify * | # * This program is free software; you can redistribute it and/or modify * | ||
| Line 25: | Line 25: | ||
# //***************************************************************************** | # //***************************************************************************** | ||
| − | + | # Data | |
| + | a=[0.,1.,2.,3.,4.] | ||
u={} | u={} | ||
| − | + | # Main loop | |
for type in ['STD','QM6']: | for type in ['STD','QM6']: | ||
u[type]=[] | u[type]=[] | ||
disp=[] | disp=[] | ||
| − | for | + | for e in a: |
# Domain | # Domain | ||
domain.planeStress(1.0) | domain.planeStress(1.0) | ||
| Line 38: | Line 39: | ||
# Nodes | # Nodes | ||
node.add( 1, 0., 0.) | node.add( 1, 0., 0.) | ||
| − | node.add( 2, 5.+ | + | node.add( 2, 5.+e, 0.) |
node.add( 3, 10., 0.) | node.add( 3, 10., 0.) | ||
node.add( 4, 0., 2.) | node.add( 4, 0., 2.) | ||
| − | node.add( 5, 5.- | + | node.add( 5, 5.-e, 2.) |
node.add( 6, 10., 2.) | node.add( 6, 10., 2.) | ||
| − | # | + | # Elements |
if type=='STD': | if type=='STD': | ||
element.quad4d( 1, 1, 2, 5, 4, 1) | element.quad4d( 1, 1, 2, 5, 4, 1) | ||
| Line 61: | Line 62: | ||
analysis.static() | analysis.static() | ||
analysis.run(1,1) | analysis.run(1,1) | ||
| − | # Get displacement | + | # Get tip displacement |
disp.append(node.data(6)['disp'][1]) | disp.append(node.data(6)['disp'][1]) | ||
domain.clear() | domain.clear() | ||
| Line 79: | Line 80: | ||
for i in range(len(a)): | for i in range(len(a)): | ||
print '|% 5.2f| % 8.4f| % 8.4f|'%(a[i],u['STD'][i],u['QM6'][i]) | print '|% 5.2f| % 8.4f| % 8.4f|'%(a[i],u['STD'][i],u['QM6'][i]) | ||
| − | print "======= | + | print "===========================" |
| − | # Plot results | + | # Plot results using matplotlib |
from pylab import * | from pylab import * | ||
plot(a,u['STD'],'b-o',linewidth=1.0,label="QM6") | plot(a,u['STD'],'b-o',linewidth=1.0,label="QM6") | ||
| Line 93: | Line 94: | ||
show() | show() | ||
</pre> | </pre> | ||
| + | [[Category:Examples]] | ||
Latest revision as of 17:16, 9 January 2008
# /******************************************************************************
# * nemesis. an experimental finite element code. *
# * Copyright (C) 2004-2008 F.E.Karaoulanis [http://www.nemesis-project.org] *
# * *
# * This program is free software; you can redistribute it and/or modify *
# * it under the terms of the GNU General Public License version 3, as *
# * published by the Free Software Foundation. *
# * *
# * This program is distributed in the hope that it will be useful, *
# * but WITHOUT ANY WARRANTY; without even the implied warranty of *
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
# * GNU General Public License for more details. *
# * *
# * You should have received a copy of the GNU General Public License *
# * along with this program. If not, see <http://www.gnu.org/licenses/>. *
# ******************************************************************************/
# //*****************************************************************************
# // $LastChangedDate$
# // $LastChangedRevision$
# // $LastChangedBy$
# // $HeadURL$
# // Author(s): F.E. Karaoulanis (fkar@nemesis-project.org)
# //*****************************************************************************
# Data
a=[0.,1.,2.,3.,4.]
u={}
# Main loop
for type in ['STD','QM6']:
u[type]=[]
disp=[]
for e in a:
# Domain
domain.planeStress(1.0)
# Material
material.elastic( 1, 1., 0.)
# Nodes
node.add( 1, 0., 0.)
node.add( 2, 5.+e, 0.)
node.add( 3, 10., 0.)
node.add( 4, 0., 2.)
node.add( 5, 5.-e, 2.)
node.add( 6, 10., 2.)
# Elements
if type=='STD':
element.quad4d( 1, 1, 2, 5, 4, 1)
element.quad4d( 2, 2, 3, 6, 5, 1)
else:
element.quad4i( 1, 1, 2, 5, 4, 1)
element.quad4i( 2, 2, 3, 6, 5, 1)
# Constraints
node.fix(1,1)
node.fix(1,2)
node.fix(4,1)
# LoadCase
lc.define(1)
load.node(3,1,+1.)
load.node(6,1,-1.)
# Analysis
analysis.static()
analysis.run(1,1)
# Get tip displacement
disp.append(node.data(6)['disp'][1])
domain.clear()
# Store displacements
u[type][:]=disp[:]
# Normalize results
u0=u['QM6'][0]
for type in ['STD','QM6']:
for i in range(len(u[type])):
u[type][i]=u[type][i]/u0
# Print results
print "==========================="
print "|alpha |STD |QM6 |"
print "+------+---------+--------+"
for i in range(len(a)):
print '|% 5.2f| % 8.4f| % 8.4f|'%(a[i],u['STD'][i],u['QM6'][i])
print "==========================="
# Plot results using matplotlib
from pylab import *
plot(a,u['STD'],'b-o',linewidth=1.0,label="QM6")
plot(a,u['QM6'],'k-o',linewidth=1.0,label="Standard")
xlabel('alpha')
ylabel('u/u0')
title('Mesh distortion test.')
ylim([0.,1.])
grid(True)
legend(loc='upper right')
show()
