Note
Go to the end to download the full example code
Patch test for triangular plate under in-plane loading
The patch test is an empirical minimum test which every finite element has to pass to ensure convergence with mesh refinement.
It consists of a problem for which a known homogeneous solution exists. For plates, we commonly use a rectangular plate subject to homogeneous edge loading, e.g., constant tension in the x-direction, or constant shear, etc.
The mesh must contain distorted elements and at least one element not attached to any node on the boundary.
Author: Peter Mackenzie-Helnwein
import numpy as np
from femedu.examples import Example
from femedu.domain import System, Node
from femedu.solver import NewtonRaphsonSolver
from femedu.elements.linear import Triangle
from femedu.materials import PlaneStress
class ExamplePlate03(Example):
def problem(self):
# ========== setting mesh parameters ==============
N = 8 # number of elements in the mesh
Lx = 100.0 # length of plate in the x-direction
Ly = 80.0 # length of plate in the y-direction
# ========== setting material parameters ==============
params = dict(
E = 20000., # Young's modulus
nu = 0.250, # Poisson's ratio
t = 1.00 # thickness of the plate
)
# ========== setting load parameters ==============
px = 10.0 # uniform load normal to x=const
py = 0.0 # uniform load normal to y=const
pxy = 0.0 # uniform shear load on x=const and y=const
# ========== setting analysis parameters ==============
target_load_level = 1.00 # reference load
max_steps = 2 # number of load steps: 2 -> [0.0, 1.0]
# define a list of target load levels
load_levels = np.linspace(0, target_load_level, max_steps)
#
# ==== Build the system model ====
#
model = System()
model.setSolver(NewtonRaphsonSolver())
# create nodes
nodes = (
Node(0.0*Lx, 0.0*Ly), # nd 0
Node(0.2*Lx, 0.0*Ly), # nd 1
Node(0.5*Lx, 0.0*Ly), # nd 2
Node(0.7*Lx, 0.0*Ly), # nd 3
Node(1.0*Lx, 0.0*Ly), # nd 4
#
Node(0.0*Lx, 0.2*Ly), # nd 5
Node(0.15*Lx,0.3*Ly), # nd 6
Node(0.5*Lx, 0.2*Ly), # nd 7
Node(0.8*Lx, 0.3*Ly), # nd 8
Node(1.0*Lx, 0.2*Ly), # nd 9
#
Node(0.0*Lx, 0.6*Ly), # nd 10
Node(0.2*Lx, 0.5*Ly), # nd 11
Node(0.7*Lx, 0.7*Ly), # nd 12
Node(0.9*Lx, 0.6*Ly), # nd 13
Node(1.0*Lx, 0.8*Ly), # nd 14
#
Node(0.0*Lx, 1.0*Ly), # nd 15
Node(0.3*Lx, 1.0*Ly), # nd 16
Node(0.55*Lx,1.0*Ly), # nd 17
Node(0.8*Lx, 1.0*Ly), # nd 18
Node(1.0*Lx, 1.0*Ly), # nd 19
)
elements = (
Triangle(nodes[0],nodes[1],nodes[5],PlaneStress(params)), # elem 0
Triangle(nodes[1],nodes[2],nodes[6],PlaneStress(params)), # elem 1
Triangle(nodes[2],nodes[3],nodes[7],PlaneStress(params)), # elem 2
Triangle(nodes[3],nodes[4],nodes[8],PlaneStress(params)), # elem 3
#
Triangle(nodes[6],nodes[5],nodes[1],PlaneStress(params)), # elem 4
Triangle(nodes[7],nodes[6],nodes[2],PlaneStress(params)), # elem 5
Triangle(nodes[8],nodes[7],nodes[3],PlaneStress(params)), # elem 6
Triangle(nodes[9],nodes[8],nodes[4],PlaneStress(params)), # elem 7
#
Triangle(nodes[5],nodes[6],nodes[10],PlaneStress(params)), # elem 8
Triangle(nodes[6],nodes[7],nodes[11],PlaneStress(params)), # elem 9
Triangle(nodes[7],nodes[8],nodes[12],PlaneStress(params)), # elem 10
Triangle(nodes[8],nodes[9],nodes[13],PlaneStress(params)), # elem 11
#
Triangle(nodes[11],nodes[10],nodes[6],PlaneStress(params)), # elem 12
Triangle(nodes[12],nodes[11],nodes[7],PlaneStress(params)), # elem 13
Triangle(nodes[13],nodes[12],nodes[8],PlaneStress(params)), # elem 14
Triangle(nodes[14],nodes[13],nodes[9],PlaneStress(params)), # elem 15
#
Triangle(nodes[10],nodes[11],nodes[15],PlaneStress(params)), # elem 16
Triangle(nodes[11],nodes[12],nodes[16],PlaneStress(params)), # elem 17
Triangle(nodes[12],nodes[13],nodes[17],PlaneStress(params)), # elem 18
Triangle(nodes[13],nodes[14],nodes[18],PlaneStress(params)), # elem 19
#
Triangle(nodes[16],nodes[15],nodes[11],PlaneStress(params)), # elem 20
Triangle(nodes[17],nodes[16],nodes[12],PlaneStress(params)), # elem 21
Triangle(nodes[18],nodes[17],nodes[13],PlaneStress(params)), # elem 22
Triangle(nodes[19],nodes[18],nodes[14],PlaneStress(params)), # elem 23
)
model.addNode(*nodes)
model.addElement(*elements)
# define support(s)
fix_x = (0,)
fix_y = (0,4)
for idx in fix_x:
nodes[idx].fixDOF('ux') # horizontal support left end
for idx in fix_y:
nodes[idx].fixDOF('uy') # vertical support right end
# ==== complete the reference load ====
# surface loads on the left side
elements[ 0].setSurfaceLoad(2,px)
elements[ 8].setSurfaceLoad(2,px)
elements[16].setSurfaceLoad(2,px)
# surface loads on the right side
elements[ 7].setSurfaceLoad(2,px)
elements[15].setSurfaceLoad(2,px)
elements[23].setSurfaceLoad(2,px)
# these are only nodal forces as part of the reference load
# .. load only the upper node
#print(model)
model.report()
model.plot(factor=0., title="undeformed system", filename="plate03_undeformed.png", show_bc=1)
model.setLoadFactor(10.0)
model.solve()
model.solver.showKt(filename="plate03_spy_Kt.png")
model.report()
model.plot(factor=100., filename="plate03_deformed.png")
Run the example by creating an instance of the problem and executing it by calling Example.run()
if __name__ == "__main__":
ex = ExamplePlate03()
ex.run()
System Analysis Report
=======================
Nodes:
---------------------
Node_646:
x: [0. 0.]
fix: ['ux', 'uy']
u: None
Node_647:
x: [20. 0.]
u: None
Node_648:
x: [50. 0.]
u: None
Node_649:
x: [70. 0.]
u: None
Node_650:
x: [100. 0.]
fix: ['uy']
u: None
Node_651:
x: [ 0. 16.]
u: None
Node_652:
x: [15. 24.]
u: None
Node_653:
x: [50. 16.]
u: None
Node_654:
x: [80. 24.]
u: None
Node_655:
x: [100. 16.]
u: None
Node_656:
x: [ 0. 48.]
u: None
Node_657:
x: [20. 40.]
u: None
Node_658:
x: [70. 56.]
u: None
Node_659:
x: [90. 48.]
u: None
Node_660:
x: [100. 64.]
u: None
Node_661:
x: [ 0. 80.]
u: None
Node_662:
x: [30. 80.]
u: None
Node_663:
x: [55. 80.]
u: None
Node_664:
x: [80. 80.]
u: None
Node_665:
x: [100. 80.]
u: None
Elements:
---------------------
Triangle_921: nodes ( Node_646 Node_647 Node_651 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_922: nodes ( Node_647 Node_648 Node_652 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_923: nodes ( Node_648 Node_649 Node_653 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_924: nodes ( Node_649 Node_650 Node_654 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_925: nodes ( Node_652 Node_651 Node_647 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_926: nodes ( Node_653 Node_652 Node_648 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_927: nodes ( Node_654 Node_653 Node_649 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_928: nodes ( Node_655 Node_654 Node_650 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_929: nodes ( Node_651 Node_652 Node_656 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_930: nodes ( Node_652 Node_653 Node_657 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_931: nodes ( Node_653 Node_654 Node_658 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_932: nodes ( Node_654 Node_655 Node_659 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_933: nodes ( Node_657 Node_656 Node_652 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_934: nodes ( Node_658 Node_657 Node_653 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_935: nodes ( Node_659 Node_658 Node_654 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_936: nodes ( Node_660 Node_659 Node_655 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_937: nodes ( Node_656 Node_657 Node_661 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_938: nodes ( Node_657 Node_658 Node_662 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_939: nodes ( Node_658 Node_659 Node_663 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_940: nodes ( Node_659 Node_660 Node_664 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_941: nodes ( Node_662 Node_661 Node_657 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_942: nodes ( Node_663 Node_662 Node_658 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_943: nodes ( Node_664 Node_663 Node_659 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
Triangle_944: nodes ( Node_665 Node_664 Node_660 )
material: PlaneStress
strain: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=-0.000e+00
stress: xx=0.000e+00 yy=0.000e+00 xy=0.000e+00 zz=0.000e+00
+
System Analysis Report
=======================
Nodes:
---------------------
Node_646:
x: [0. 0.]
fix: ['ux', 'uy']
u: [0. 0.]
Node_647:
x: [20. 0.]
u: [ 1.00000000e-01 -1.42162675e-16]
Node_648:
x: [50. 0.]
u: [2.50000000e-01 5.58828228e-16]
Node_649:
x: [70. 0.]
u: [ 3.50000000e-01 -4.75859138e-15]
Node_650:
x: [100. 0.]
fix: ['uy']
u: [0.5 0. ]
Node_651:
x: [ 0. 16.]
u: [-1.02467033e-15 -2.00000000e-02]
Node_652:
x: [15. 24.]
u: [ 0.075 -0.03 ]
Node_653:
x: [50. 16.]
u: [ 0.25 -0.02]
Node_654:
x: [80. 24.]
u: [ 0.4 -0.03]
Node_655:
x: [100. 16.]
u: [ 0.5 -0.02]
Node_656:
x: [ 0. 48.]
u: [ 1.37573745e-15 -6.00000000e-02]
Node_657:
x: [20. 40.]
u: [ 0.1 -0.05]
Node_658:
x: [70. 56.]
u: [ 0.35 -0.07]
Node_659:
x: [90. 48.]
u: [ 0.45 -0.06]
Node_660:
x: [100. 64.]
u: [ 0.5 -0.08]
Node_661:
x: [ 0. 80.]
u: [-2.5888467e-15 -1.0000000e-01]
Node_662:
x: [30. 80.]
u: [ 0.15 -0.1 ]
Node_663:
x: [55. 80.]
u: [ 0.275 -0.1 ]
Node_664:
x: [80. 80.]
u: [ 0.4 -0.1]
Node_665:
x: [100. 80.]
u: [ 0.5 -0.1]
Elements:
---------------------
Triangle_921: nodes ( Node_646 Node_647 Node_651 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-7.115e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=1.187e-12 xy=-5.692e-13 zz=0.000e+00
Triangle_922: nodes ( Node_647 Node_648 Node_652 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-6.687e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=-2.370e-12 xy=-5.349e-13 zz=0.000e+00
Triangle_923: nodes ( Node_648 Node_649 Node_653 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-2.659e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-1.183e-12 xy=-2.127e-12 zz=0.000e+00
Triangle_924: nodes ( Node_649 Node_650 Node_654 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-1.050e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-1.183e-12 xy=-8.400e-13 zz=0.000e+00
Triangle_925: nodes ( Node_652 Node_651 Node_647 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-1.388e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=1.187e-12 xy=-1.110e-12 zz=0.000e+00
Triangle_926: nodes ( Node_653 Node_652 Node_648 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=5.551e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=-1.183e-12 xy=4.441e-13 zz=0.000e+00
Triangle_927: nodes ( Node_654 Node_653 Node_649 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=1.665e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-1.183e-12 xy=1.332e-12 zz=0.000e+00
Triangle_928: nodes ( Node_655 Node_654 Node_650 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-5.551e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=1.183e-12 xy=-4.441e-13 zz=0.000e+00
Triangle_929: nodes ( Node_651 Node_652 Node_656 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-3.601e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=-1.183e-12 xy=-2.881e-13 zz=0.000e+00
Triangle_930: nodes ( Node_652 Node_653 Node_657 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-2.776e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=-2.370e-12 xy=-2.220e-13 zz=0.000e+00
Triangle_931: nodes ( Node_653 Node_654 Node_658 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=1.110e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=4.736e-12 xy=8.882e-13 zz=0.000e+00
Triangle_932: nodes ( Node_654 Node_655 Node_659 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=3.886e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=2.370e-12 xy=3.109e-12 zz=0.000e+00
Triangle_933: nodes ( Node_657 Node_656 Node_652 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-3.331e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-4.739e-12 xy=-2.665e-12 zz=0.000e+00
Triangle_934: nodes ( Node_658 Node_657 Node_653 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-5.551e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=1.187e-12 xy=-4.441e-13 zz=0.000e+00
Triangle_935: nodes ( Node_659 Node_658 Node_654 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=2.776e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=2.370e-12 xy=2.220e-12 zz=0.000e+00
Triangle_936: nodes ( Node_660 Node_659 Node_655 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=5.836e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=1.183e-12 xy=4.668e-12 zz=0.000e+00
Triangle_937: nodes ( Node_656 Node_657 Node_661 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-2.904e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-3.553e-12 xy=-2.323e-12 zz=0.000e+00
Triangle_938: nodes ( Node_657 Node_658 Node_662 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=2.220e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-3.553e-12 xy=1.776e-12 zz=0.000e+00
Triangle_939: nodes ( Node_658 Node_659 Node_663 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-2.220e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-2.366e-12 xy=-1.776e-12 zz=0.000e+00
Triangle_940: nodes ( Node_659 Node_660 Node_664 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=1.943e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=1.187e-12 xy=1.554e-12 zz=0.000e+00
Triangle_941: nodes ( Node_662 Node_661 Node_657 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=2.599e-17 zz=-9.375e-04
stress: xx=1.000e+02 yy=-1.183e-12 xy=2.079e-13 zz=0.000e+00
Triangle_942: nodes ( Node_663 Node_662 Node_658 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-1.135e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-9.475e-12 xy=-9.078e-13 zz=0.000e+00
Triangle_943: nodes ( Node_664 Node_663 Node_659 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=-3.331e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=-3.553e-12 xy=-2.665e-12 zz=0.000e+00
Triangle_944: nodes ( Node_665 Node_664 Node_660 )
material: PlaneStress
strain: xx=5.000e-03 yy=-1.250e-03 xy=7.105e-16 zz=-9.375e-04
stress: xx=1.000e+02 yy=8.288e-12 xy=5.684e-12 zz=0.000e+00
/usr/local/lib/python3.11/site-packages/matplotlib/quiver.py:632: RuntimeWarning: Mean of empty slice.
amean = a.mean()
/usr/local/lib/python3.11/site-packages/numpy/core/_methods.py:129: RuntimeWarning: invalid value encountered in scalar divide
ret = ret.dtype.type(ret / rcount)
Total running time of the script: (0 minutes 0.548 seconds)