1# =============================================================================
2# PROJECT CHRONO - http:#projectchrono.org
3#
4# Copyright (c) 2014 projectchrono.org
5# All rights reserved.
6#
7# Use of this source code is governed by a BSD-style license that can be found
8# in the LICENSE file at the top level of the distribution and at
9# http:#projectchrono.org/license-chrono.txt.
10#
11# =============================================================================
12# Authors: Simone Benatti
13# =============================================================================
14#
15# FEA demo on applying loads to beams
16#
17# =============================================================================
18
19import pychrono as chrono
20import pychrono.fea as fea
21import pychrono.irrlicht as chronoirr
22import os
23import copy
24
25out_dir = chrono.GetChronoOutputPath() + "FEA_LOADS"  # Output directory
26
27
28print("Copyright (c) 2017 projectchrono.org ")
29
30# Create (if needed) output directory
31if not os.path.isdir(out_dir):
32    os.mkdir(out_dir)
33
34
35# Create the physical system
36my_system = chrono.ChSystemSMC()
37
38# Create the Irrlicht visualization
39
40application = chronoirr.ChIrrApp(my_system, "Loads on beams", chronoirr.dimension2du(800, 600))
41application.AddTypicalLogo()
42application.AddTypicalSky()
43application.AddTypicalLights()
44application.AddTypicalCamera(chronoirr.vector3df(0.5, 0.0, -3.0), chronoirr.vector3df(0.5, 0.0, 0.0))
45
46# Create a mesh
47mesh = fea.ChMesh()
48my_system.Add(mesh)
49
50# Create some nodes (with default mass 0)
51nodeA = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(0, 0, 0)))
52nodeB = fea.ChNodeFEAxyzrot(chrono.ChFrameD(chrono.ChVectorD(2, 0, 0)))
53nodeA.SetMass(0.0)
54nodeB.SetMass(0.0)
55mesh.AddNode(nodeA)
56mesh.AddNode(nodeB)
57
58# Create beam section & material
59msection = fea.ChBeamSectionEulerAdvanced()
60beam_wy = 0.1
61beam_wz = 0.2
62msection.SetAsRectangularSection(beam_wy, beam_wz)
63msection.SetYoungModulus(0.01e9)
64msection.SetGshearModulus(0.01e9 * 0.3)
65msection.SetBeamRaleyghDamping(0.200)
66msection.SetDensity(1500)
67
68# Create an Eulero-Bernoulli beam with a single element
69elementA = fea.ChElementBeamEuler()
70elementA.SetNodes(nodeA, nodeB)
71elementA.SetSection(msection)
72mesh.AddElement(elementA)
73
74# Create the ground body
75ground = chrono.ChBody()
76ground.SetBodyFixed(True)
77my_system.Add(ground)
78
79# Create a constraat the end of the beam
80constrA = chrono.ChLinkMateGeneric()
81constrA.Initialize(nodeA, ground, False, nodeA.Frame(), nodeA.Frame())
82my_system.Add(constrA)
83constrA.SetConstrainedCoords(True, True, True,   # x, y, z
84                              True, True, True)  # Rx, Ry, Rz
85
86# -----------------------------------------------------------------
87# Apply loads
88
89# Create the load container
90loadcontainer = chrono.ChLoadContainer()
91my_system.Add(loadcontainer)
92
93
94#define LOAD_5
95
96
97print("Applied loads: \n")
98
99
100
101# Example 6:
102# As before, create a custom load with stiff force, acting on a single node, but
103# this time we inherit directly from ChLoadCustom, i.e. a load that does not require ChLoader features.
104# This is mostly used in case one does not need the automatic surface/volume quadrature of ChLoader.
105# As a stiff load, this will automatically generate a jacobian (tangent stiffness matrix K)
106# that will be used in statics, implicit integrators, etc.
107
108print("   Custom load with stiff force, acting on a single node (VER 2).")
109
110nodeD = fea.ChNodeFEAxyz(chrono.ChVectorD(2, 10, 3))
111mesh.AddNode(nodeD)
112
113class MyLoadCustom(chrono.ChLoadCustom):
114    def __init__(self, mloadable):
115        chrono.ChLoadCustom.__init__(self, mloadable)
116    #/ "Virtual" copy constructor (covariant return type).
117    def Clone(self):
118        newinst = copy.deepcopy(self)
119        return  newinst
120    # Compute Q=Q(x,v)
121    # This is the function that you have to implement. It should return the generalized Q load
122    # (i.e.the force in generalized lagrangian coordinates).
123    # For ChNodeFEAxyz, Q loads are expected as 3-rows vectors, containing absolute force x,y,z.
124    # As this is a stiff force field, dependency from state_x and state_y must be considered.
125    def ComputeQ(self,state_x,      #/< state position to evaluate Q
126                 state_w):     #/< state speed to evaluate Q
127        if not state_x==None and not state_w==None :
128            node_pos = chrono.ChVectorD(state_x[0], state_x[1], state_x[2])
129            node_vel = chrono.ChVectorD(state_w[0], state_w[1], state_w[2])
130        else:
131            mynode = fea.CastToChNodeFEAxyz( fea.CastToChNodeFEAbase( chrono.CastToChNodeBase(self.loadable) ))
132            node_pos = mynode.GetPos()
133            node_vel = mynode.GetPos_dt()
134        # Just implement a simple force+spring+damper in xy plane,
135        # for spring&damper connected to absolute reference:
136        Kx = 100
137        Ky = 400
138        Dx = 0.6
139        Dy = 0.9
140        x_offset = 2
141        y_offset = 5
142        x_force = 50
143        y_force = 0
144        # Store the computed generalized forces in this.load_Q, same x,y,z order as in state_w
145        self.load_Q[0] = x_force - Kx * (node_pos.x - x_offset) - Dx * node_vel.x
146        self.load_Q[1] = y_force - Ky * (node_pos.y - y_offset) - Dy * node_vel.y
147        self.load_Q[2] = 0
148    # Remember to set this as stiff, to enable the jacobians
149    def IsStiff(self) :
150        return True
151
152# Instance load object, applying to a node, as in previous example, and add to container:
153mloadcustom = MyLoadCustom(nodeD)
154loadcontainer.Add(mloadcustom)
155
156
157
158# Stiffness and damping matrices
159K_matrix = chrono.ChMatrixDynamicD(6, 6)
160R_matrix = chrono.ChMatrixDynamicD(6, 6)
161
162for i in range(6):
163    K_matrix[i, i] = 1e5
164
165for i in range(6):
166    for j in range(6):
167        R_matrix[i, j] = 1e-3
168
169ch_bushing = fea.ChLoadXYZROTnodeBodyBushingGeneric(
170    nodeB,
171    ground,
172    nodeB,
173    K_matrix,
174    R_matrix)
175
176loadcontainer.Add(ch_bushing)
177
178# -----------------------------------------------------------------
179# Set visualization of the FEM mesh.
180
181mvisualizebeamA = fea.ChVisualizationFEAmesh(mesh)
182mvisualizebeamA.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_ELEM_BEAM_MZ)
183mvisualizebeamA.SetColorscaleMinMax(-400, 200)
184mvisualizebeamA.SetSmoothFaces(True)
185mvisualizebeamA.SetWireframe(False)
186mesh.AddAsset(mvisualizebeamA)
187
188mvisualizebeamC = fea.ChVisualizationFEAmesh(mesh)
189mvisualizebeamC.SetFEMglyphType(fea.ChVisualizationFEAmesh.E_GLYPH_NODE_CSYS)
190mvisualizebeamC.SetFEMdataType(fea.ChVisualizationFEAmesh.E_PLOT_NONE)
191mvisualizebeamC.SetSymbolsThickness(0.006)
192mvisualizebeamC.SetSymbolsScale(0.01)
193mvisualizebeamC.SetZbufferHide(False)
194mesh.AddAsset(mvisualizebeamC)
195
196application.AssetBindAll()
197application.AssetUpdateAll()
198
199# -----------------------------------------------------------------
200
201# Setup a MINRES solver. For FEA one cannot use the default PSOR type solver.
202solver = chrono.ChSolverMINRES()
203my_system.SetSolver(solver)
204solver.SetMaxIterations(200)
205solver.SetTolerance(1e-15)
206solver.EnableDiagonalPreconditioner(True)
207solver.SetVerbose(False)
208
209my_system.SetSolverForceTolerance(1e-13)
210
211# Set integrator
212ts = chrono.ChTimestepperEulerImplicitLinearized(my_system)
213my_system.SetTimestepper(ts)
214
215# Simulation loop
216
217application.SetTimestep(1e-3)
218
219while (application.GetDevice().run()) :
220    application.BeginScene()
221    application.DrawAll()
222    application.DoStep()
223    if not application.GetPaused() :
224        time = my_system.GetChTime()
225        posB = nodeB.GetPos()
226    application.EndScene()
227
228