User’s Guide¶
This guide introduces the components and commands that define the configuration file required to execute an analysis using MYFEMPY. Depending on the problem setup, adjustments may be necessary; please refer to the help or API documentation for details. Additional Tutorials are provided at the end of this guide.
Presentation¶
The myfempy project was developed as a finite element solver for multi-physics simulations, intended for academic and research use. Originating from the author’s doctoral thesis, the project has matured into a robust computational framework. Its primary goal is to provide an accessible yet powerful environment for simulation setup and execution, minimizing the need for advanced prior knowledge of numerical methods while maintaining flexibility for expert users.
At its core, myfempy implements a finite element processor that drives the numerical analysis. The workflow is organized into three principal stages:
Pre-processing
- Configuration is performed through API commands.
- Users specify material properties, element types, and boundary conditions.
- Mesh generation is managed via Gmsh, which supplies both geometry and discretization data.
Processing
- The solver core executes finite element computations according to the defined physics and constraints.
- Multi-physics coupling is supported, enabling integration of multiple physical domains within a single simulation.
Post-processing
- Results are exported in formats compatible with visualization tools.
- Output data is structured as VTK files, which can be analyzed and visualized using ParaView.
This architecture emphasizes modularity, extensibility, and interoperability with established open-source tools, positioning myfempy as a versatile platform for computational mechanics and multi-physics research. A typical analysis workflow involves submitting requests (Python dictionaries) through a specific script pattern. The user imports the desired solver, configures the problem by defining materials, finite elements, and geometry, and generates a mesh via Gmsh. Loads and boundary conditions are applied, and the solver computes the solution. Results are exported and visualized in ParaView.
graph LR
U[User] --> MyfempyAPI
MyfempyAPI --> A[Gmsh]
A[Gmsh] --> B[MyfempyCore]
B --> C[Paraview]
style A fill:#006400, stroke:#333, stroke-width:2px
style B fill:#006400, stroke:#333, stroke-width:2px
style C fill:#00008B, stroke:#333, stroke-width:2pxTo ensure organized execution, myfempy employs a hierarchical code structure that integrates user requests into a processing pipeline. Input data is validated through controllers in the I/O interface. If commands and data are consistent, they are passed to the solver core, where pre-processing, solution, and post-processing occur. The processed data is then returned through the I/O filters and exported as visualization-ready files (VTK).
graph
U[InputData] --> B{IO}
B -->|Yes| C[PreProcess]
C --> R[CoreSolver]
R --> D[PostProcess]
D --> B;
B -----> O[OutputData]
B ----->|No| E[ReturnERROR]
style C fill:#006400, stroke:#333, stroke-width:2px
style R fill:#006400, stroke:#333, stroke-width:2px
style D fill:#006400, stroke:#333, stroke-width:2px
style O fill:#00008B, stroke:#333, stroke-width:2px
style E fill:#8B0000,stroke:#333,stroke-width:2px;How it Works¶
myfempy is implemented using an object-oriented paradigm, structured around the Bridge design pattern. The system is composed of multiple classes, with user interaction facilitated through APIs that orchestrate the execution of the solver.
The project adopts a black-box architecture: users provide input data, the solver executes the numerical routines, and the system returns output data along with logs for validation.
graph LR
UserInput --> BLACK-BOX
BLACK-BOX --> OutputDataThe main API provides direct access to the core components of the project, enabling construction of the Model, definition of Physics, and execution of the Solution.
graph LR
User --> API
API --> Model
API --> Physics
API --> Solve
API --> Results
Model --> Element
Model --> Shape
Model --> Mesh
Model --> Material
Model --> Geometry
Physics --> Loads
Physics --> BoundCond
Physics --> Coupling
Results --> PostProcess
Solve --> Assembler
Solve --> SolverPre-Process¶
API and Solver¶
# ==============================================================================
# import newAnalysis to set a new problem
# ==============================================================================
from myfempy import newAnalysis
# ==============================================================================
# import a solver to running and analysis the problem set, e.g. SteadyStateLinear
# ==============================================================================
from myfempy import Solver
# ==============================================================================
# now, fea is the API to running the Solver Analysis
# ==============================================================================
fea = newAnalysis(Solver)
Solvers¶
# ==============================================================================
# currently available solvers for import into the myfempy library
# ==============================================================================
# >see: Table 3 - Solvers List to more informations
# solver options
"SteadyStateLinear" # Steady State Linear Solver Class
"SteadyStateLinearIterative" # Steady State Linear Iterative Solver Class
"DynamicEigenLinear" # Dynamic Eigen (modal problem) Linear Solver Class
"DynamicHarmonicResponseLinear" # Dynamic Harmonic Response Forced System Steady State Linear Solver Class
# [adv]
"StaticLinearCyclicSymmPlane" # Static Linear Cyclic Symmetry Plane Solver Class
"PhononicCrystalInPlane" # Phononic Crystal In-Plane Solver Class
"HomogenPlane" # Homogenization Plane Solver Class
Model¶
# ==============================================================================
# modeldata is a Python dictionary that contains the commands for model set
# ==============================================================================
modeldata = dict()
Material¶
# ==============================================================================
# config. material
# ==============================================================================
modeldata["MATERIAL"] = dict()
keys¶
"MAT":str() # material set
# >see: Table 4 - Material List to more informations
# options
'lumped' # lumped material
'uniaxialstress' # axial{rod, beams...} behavior material
'planestress' # plane stress behavior
'planestrain' # plane strain behavior
'solidelastic' # solid behavior material
'heatplane' # heat behavior material
'heatsolid' # heat behavior material
# [dev]
'fluid'
"TYPE":str() # material behavior
# options
'isotropic' # isotropic stress/strain material
'usermaterial' # user config. material
"PROPMAT":list[mat_set_1:dict(),..., mat_set_n:dict()] # material properties
# options
"NAME": str() # material name def
# isotropic solid parameters
# >see: Appendix: Table 5 - Consistent Units to more informations
"EXX": float() # elasticity modulus in x direction
"VXY": float() # poisson's ratio in x direction
"GXY": float() # shear modulus in x direction
"RHO": float() # density of material
# anisotropic solid parameters
"EYY": float()
"VYZ": float()
"GYZ": float()
"EZZ": float()
"VZX": float()
"GZX": float()
# heat parameters
"KXX": float() # thermal conductivity in x direction
"KYY": float() # thermal conductivity in y direction
"KZZ": float() # thermal conductivity in z direction
"CTE": float() # coefficients thermal expansion
# [dev] fluid parameters
"VIS": float() # viscosity
"RHO": float() # density of fluid
# lumped model
"STIF": float() # stiffness lumped
"DAMP": float() # damping lumped
# [adv]
"CLASS": class() # user new class defined
#-----------------------
# example:
#-----------------------
# from myfempy import PlaneStress
# ===============================================================================
# SET NEW MATERIAL <ORTHOTROPIC ELASTIC>
# ===============================================================================
# class UserNewMaterial(PlaneStress):
# def __init__(self):
# super().__init__()
# def getMaterialSet():
# matset = {
# "mat": "planestress",
# "type": "orthotropic",
# }
# return matset
# def getElasticTensor(tabmat, inci, element_number, a=None, b=None):
# # # material elasticity
# EXX = tabmat[int(inci[element_number, 2]) - 1]["EXX"]
# EYY = tabmat[int(inci[element_number, 2]) - 1]["EYY"]
# # material poisson ratio
# VXY = tabmat[int(inci[element_number, 2]) - 1]["VXY"]
# VYZ = tabmat[int(inci[element_number, 2]) - 1]["VYZ"]
# # EXX = 45E3 # N/mm^2 --> MPa
# # EYY = 12E3 # N/mm^2 --> MPa
# # VXY = 0.23
# # VYX = 0.66
# S00 = EXX/(1-VXY*VYZ)
# S01 = (VYZ*EXX)/(1-VXY*VYZ)
# S10 = (VXY*EYY)/(1-VXY*VYZ)
# S11 = EYY/(1-VXY*VYZ)
# S22 = (EXX*EYY)/(EXX+EYY+2.0*EYY*VXY)
# D = np.zeros((3, 3), dtype=np.float64)
# D[0, 0] = S00
# D[0, 1] = S01
# D[1, 0] = S10
# D[1, 1] = S11
# D[2, 2] = S22
# return D
# def getFailureCriteria(stress):
# stress
# return super().getFailureCriteria()
# ===============================================================================
Geometry¶
# ==============================================================================
# config. geometry
# ==============================================================================
modeldata["GEOMETRY"] = dict()
keys¶
"GEO":str() # geometry type
# options
'solid' # solid 3D geo.
'thickness' # planes geometries, e.g. plane stress, plates ...
'frame' # rods, beams and frames geometries
"SECTION":str() # geometry beam cross section set
# options
# >see: Appendix: Cross Section Dimensions to more informations
'rectangle'
'rectangle_tube'
'circle'
'circle_tube'
'isection'
'tsection'
'csection'
'lsection'
'userdefined'
"PROPGEO":list[geo_set_1:dict(),..., geo_set_n:dict()] # geometry properties
# options
"NAME": str() # geometry name def
# parameters defined from the user's cross-section
"AREACS": float() # area cross section
"INERXX":float() # inercia x diretion
"INERYY":float() # inercia y diretion
"INERZZ":float() # inercia z diretion
"THICKN":float() # thickness of plane/plate
# user set section dimensions
# >see: Appendix: Cross Section Dimensions to more informations
"DIM":dict() # dimensional beam cross section
{
"b":float() # b size
"h":float() # h size
"t":float() # t size
"d":float() # d size
}
# user set Center of Gravity coord. section
"CG":dict() # center of gravity of beam cross section
{
"y_max":float(),
"y_min":float(),
"z_max":float(),
"z_min":float(),
"r_max":float()
}
Element¶
# ==============================================================================
# config. element
# ==============================================================================
modeldata["ELEMENT"] = dict()
keys¶
"TYPE": str() # finite element
# options
# >see: Table 1 - Elements List to more informations
'structbeam' # Beam Structural
'structplane' # Plane Structural
'structsolid' # Solid Structural
'heatplane' # Plane Heat
'heatsolid' # Solid Heat
# [dev]
'structplate'
'fluid2d'
'fluid3d'
'SHAPE': str()
# options
# >see: Table 2 - Mesh List to more informations
'line2' # Line 2-Node
'line3' # Line 3-Node
'tria3' # Triangular 3-Node
'tria6' # Triangular 6-Node
'quad4' # Quadrilateral 4-Node
'quad8' # Quadrilateral 8-Node
'hexa8' # Hexaedron 8-Node
'tetr4' # Tetrahedron 4-Node
'usershape' # User Defined Shape
# optional
'INTGAUSS': int() # number of points to integrations the mesh
# >see: Table 2 - Mesh List to more informations
Mesh¶
# ==============================================================================
# config. mesh
# ==============================================================================
modeldata["MESH"] = dict()
Manual Mesh Options¶
'TYPE': 'manual'
'COORD': nodes_coord_array
# [
# [node_number_1:int, coord_x_node_1:float, coord_y_node_1:float, coord_z_node_1:float]
# ...
# [node_number_N:int, coord_x_node_N:float, coord_y_node_N:float, coord_z_node_N:float]
# ]
#-----------------------
# example:
#-----------------------
# >>nodes_coord_array =
# [
# [1, 0, 0, 0]
# [2, 1, 0, 0]
# [3, 0, 1, 0]
# ]
'INCI': mesh_incidence_array
# [
# [elem_number_n:int, mat_type:int, geo_type:int, nodes_list_conec_1,...nodes_list_conec_n]
# ...
# ]
#-----------------------
# example:
#-----------------------
# >>mesh_incidence_array =
# [[1, 1, 1, 1, 2, 3, 4]]
Legacy Mesh Options¶
'TYPE': 'legacy'
'LX': float() # set a length in x diretion
'LY': float() # set a length in y diretion
'NX': int() # set a number of elements in x diretion
'NY': int() # set a number of elements in y diretion
Gmsh Mesh Options [adv]¶
'TYPE': 'gmsh',
'meshimport':dict{} # opt. to import a external gmsh mesh
# option
'object':str(object name [.msh2]) # file .msh2 only, mesh from gmsh [current version]
'cadimport':dict{} # opt. to import a cad model from any cad program
# option
'object':str(object name [.step]) # file .step/.stp only [current version]
# Options to build a self model in .geo file (from gmsh)
'filename':str() # names of the output files generated by the internal gmsh generator (.geo, .msh2)
'pointlist': list[] # poinst coord. list
# set
# [
# [coord_x_point_1:float, coord_y_point_1:float, coord_z_point_1:float]
# ...
# ]
# y
# |
# |
# (1)----x
# \
# \
# z
#-----------------------
# example:
#-----------------------
# points = [
# [0, 0, 0],
# [10, 0, 0],
# [10, 20, 0],
# [0, 20, 0]
# ]
'linelist': list[] # lines points conec., counterclockwise count
# set
# [
# [point_i_line_1:int, point_j_line_1:int]
# ...
# ]
# (i)-----{1}-----(j)
#-----------------------
# example:
#-----------------------
# lines = [
# [1, 2],
# [2, 3],
# [3, 4],
# [4, 1],
# ]
'circle':list[] # circle line, counterclockwise count
# set
# [
# [R,[CX,CY,CZ],[A0, A1]] # arc_1
# ...
# ]
# A1 ^
# | /
# | /
# | R
# | /
# |/
# (i:CX,CY,CZ)------A0
# options
R:float() # circle radius
CX:float() # point i center x coord.
CY:float() # point i center y coord.
CZ:float() # point i center z coord.
A0:str(e.g. val.='0') # angle begin [rad]
A1:str(e.g. val.='Pi/2') # angle end [rad]
#-----------------------
# example:
#-----------------------
# circle = [[30, [100, 100, 0], ['0', '2*Pi']]
# circles have priority numbering over arcs; pay attention to the plane numbering.
# see the generated model (.geo) in gmsh for correct setting
'arc': list[] # arc 3 points needed
# set
# [
# [point_i_begin:int, point_j_midle:int, , point_k_end:int]
# ...
# ]
#-----------------------
# example:
#-----------------------
# points = [
# [0, 0, 0], # ponto 1
# [200, 0, 0], # ponto 2
# [200, 40, 0], # ponto 3
# [100, 40, 0], # ponto 4
# [90, 30, 0], # ponto 5
# [0, 30, 0], # ponto 6
# [90, 40, 0], # ponto 7
# ]
# lines = [
# [1, 2], # linha 1
# [2, 3], # linha 2
# [3, 4], # linha 3
# [5, 6], # linha 4
# [6, 1], # linha 5
# ]
# arcs = [
# [5, 7, 4], # arco, adiciona uma nova linha 6
# ]
# plane = [
# [1, 2, 3, 6, 4, 5], # plane 1, sequencia das linhas
# ]
# You can add a list of lines to form a plan.
# To remove a plan within another, add the "-" symbol to the number.
# If the same number is followed by another sequence with a positive value,
# the program will remove and add a plan, creating a new material.
'planelist': list[] # planes lines conec., counterclockwise count
# set
# [
# [line_1_plane_1:int, ..., line_n_plane_1:int]
# ...
# ]
# (l)-----{3}-----(k)
# | |
# | |
# {4} [1] {2}
# | |
# | |
# (i)-----{1}-----(j)
#-----------------------
# example:
#-----------------------
# plane = [[1, 2, 3, 4]] add a new plane with the lines 1, 2, 3, 4
# plane = [[1, 2, 3, 4], [-5]] remove the line 5 ('circle') from the main plane with the lines 1, 2, 3, 4
# plane = [[1, 2, 3, 4], [-5], [5]] remove the line 5 ('circle') and add a new plane 5 ('circle') from the main plane with the lines 1, 2, 3, 4
'meshconfig':dict{} # mesh configuration inputs
# options
'mesh': str() # !!! mandatory option !!! set a type of mesh used in analysis
'sizeelement':float() # size min. of elements
'numbernodes':int() # select a number of nodes in line, only to 'line2/ line3' >see: Table 2 - Mesh List to more informations
'extrude':float() # extrude dimensional, in z diretion, from a xy plane
'meshmap':dict{} # gen. a mapped structured mesh
# option
'on': bool() # turn on(true/ false)
True
False
'edge': [ # select edge(lines) to map (only in 'on':True)
'all'
or
TAGS NUMB:list[int()] # select all edge or a specific edge
]
'numbernodes':list[int()] # select a number of nodes in each edge
#-----------------------
# example:
#-----------------------
# 'meshmap': {'on': True,
# 'edge': [[1, 2, 3], [4, 5, 6, 7]] or 'all',
# "numbernodes": [12, 8],
# [ adv ]
# Developed a robust Gmsh script that reorders nodes based on spatial coordinates (Z-Y-X)
# while preserving all physical groups and entity metadata for optimized FEM assembly.
# It is not compatible with 'meshmap'.
'reordermesh': bool() # turn on(true/ false)
True
False
# ==============================================================================
# Finally, pass the modeldata to the Model constructor API with the configured data and commands.
# ==============================================================================
fea.Model(modeldata)
Physics Setting¶
# ==============================================================================
# physicdata is a Python dictionary that contains the commands for physics set
# ==============================================================================
physicdata = dict()
# ==============================================================================
# physic set
# ==============================================================================
physicdata["PHYSIC"] = dict()
Domain¶
# ==============================================================================
# set the domain
# ==============================================================================
physicdata["PHYSIC"]["DOMAIN"]: str()
# options
'structural'
'thermal'
# [dev]
'fluid'
Loads¶
# ==============================================================================
# configuration of loads applied to physics
# ==============================================================================
physicdata["PHYSIC"]["LOAD"] = list[load_set_1:dict(),..., load_set_n:dict()]
keys¶
"TYPE":str() # type force n def.
# options
'forcenode' # force in nodes, concentrated load
'forceedge' # force in edge, distributed load
'forcesurf' # force in surface, distributed load
'forcebeam' # force in beam only opt., distributed load
# [adv]
'forcebody'
'strainzero'
# heat options
'heatfluxedge'
'heatfluxsurf'
'convectionedge'
'convectionsurf'
'heatgeneration'
"DOF":str() # dof direction of force n
# options
'fx' # force in x dir.
'fy' # force in y dir.
'fz' # force in z dir.
'pressure' # pressure app. in normal direction
'tx' # torque/moment in x dir.
'ty' # torque/moment in y dir.
'tz' # torque/moment in z dir.
'masspoint' # mass concentrated applied in node/point
'spring2ground' # spring connected node to ground/fixed end
'damper2ground' # damper connected node to ground/fixed end
# ----- OPT. WITH LOC SEEKERS
"DIR":str() # set direction >see: Axis Diretions
# options
'node' # node in mesh
'lengthx' # length line in x dir., beam only option [legacy version]
'lengthy' # length line in y dir., beam only option [legacy version]
'lengthz' # length line in z dir., beam only option [legacy version]
'edgex' # edge def in x dir. >'LOC': {'x':float(coord. x nodes), 'y':999(select all node in y dir.), 'z':float(coord. z nodes)}
'edgey' # edge def in y dir.
'edgez' # edge def in z dir.
'surfxy' # surf def in xy plane >'LOC': {'x':999, 'y': 999, 'z':float(coord. z nodes)}
'surfyz' # surf def in yz plane
'surfzx' # surf def in zx plane
"LOC":dict() # coord. node locator >see: Axis Diretions
'x':float() # x coord. node
'y':float() # y coord. node
'z':float() # z coord. node
# ----- OPT. WITH TAG SEEKERS
"DIR":str() # set direction >see: Axis Diretions
# options
'point' # point number in tag list
'edge' # edge number in tag list
'surf' # surface number in tag list
"TAG":int() # tag number of regions type, used with gmsh mesh gen, view list
# look up the corresponding node number in the Paraview, and add 1.
# For example: node_paraview = 0 --> 0 + 1 --> MESHNODE = [1]
"MESHNODE": list[int()] # select the node ([1]) or nodes list([1, 2, 3, ...]) number directly from the mesh
"VAL":list() # value list of force on steps, signal +/- is the direction
Boundary Conditions¶
# ==============================================================================
# configuration of boundary conditions applied to physics
# ==============================================================================
physicdata["PHYSIC"]["BOUNDCOND"] = list[bc_set_1:dict(),..., bc_set_n:dict()]
keys¶
"TYPE":str() # type force n def.
# options
'fixed' # fixed boundary condition u=0. More in
'displ' # displ boundary condition u!=0.
# [adv]
'cycsym'
'bloch'
# heat options
'insulated'
'temperature'
"DOF":str() # dof direction of force n
# options
'ux' # displ in x dir.
'uy' # displ in y dir.
'uz' # displ in z dir.
'rx' # rotation in x dir.
'ry' # rotation in y dir.
'rz' # rotation in z dir.
'full' # full dof
# ----- OPT. WITH LOC SEEKERS
"DIR":str() # set direction >see: Axis Diretions
# options
'node' # node in mesh
'edgex' # edge def in x dir. >'LOC': {'x':float(coord. x nodes), 'y':999(select all node in y dir.), 'z':float(coord. z nodes)}
'edgey' # edge def in y dir.
'edgez' # edge def in z dir.
'surfxy' # surf def in xy plane >'LOC': {'x':999, 'y': 999, 'z':float(coord. z nodes)}
'surfyz' # surf def in yz plane
'surfzx' # surf def in zx plane
"LOC":dict()
# options # coord. node locator >see: Axis Diretions
'x':float() # x coord. node
'y':float() # y coord. node
'z':float() # z coord. node
# ----- OPT. WITH TAG SEEKERS
"DIR":str() # set direction >see: Axis Diretions
# options
'point' # point number in tag list
'edge' # edge number in tag list
'surf' # surface number in tag list
"TAG":int() # tag number of regions type, used with gmsh mesh gen, view list
# look up the corresponding node number in the Paraview, and add 1.
# For example: node_paraview = 0 --> 0 + 1 --> MESHNODE = [1]
"MESHNODE": list[int()] # select the node ([1]) or nodes list([1, 2, 3, ...]) number directly from the mesh
"VAL":list() # value list of dislp on steps
Coupling¶
# ==============================================================================
# configuration of coupling applied to physics
# ==============================================================================
physicdata["COUPLING"] = dict()
keys¶
"TYPE":
# options
'thermalstress', # CouplSolver
# [dev]
'fsi'
# To perform a coupled analysis, the same model is used, and the post-processing
# data from the first simulation is passed as input to the coupled simulation.
"POST": [PostProcessData] # post-processing data from the previous simulation
# myfempy solver blackbox structure
# Fist Simulation [1]
# input API outputs
# Model [1] +---------------+
# BounCond + Loads [1] --> | FistSolver | --> PostProcessData [1]
# SolverSet [1] +---------------+
# Coupled Simulation [2]
# input API outputs
# Model [1] +---------------+
# BounCond + Loads [2] --> | CouplSolver | --> PostProcessData [2]
# PostProcessData [1] +---------------+
# SolverSet [2]
# ==============================================================================
# Finally, pass the physicdata to the Physic constructor API with the configured data and commands
# ==============================================================================
fea.Physic(physicdata)
Preview analysis¶
# ==============================================================================
# Preview set
# ==============================================================================
previewset = {'RENDER':
{
'filename': str(),
'show': bool()
# options
True
False
'scale': int(),
'savepng': bool(),
# options
True
False
'lines': bool(), # wireframe lines
# options
True
False
'plottags': {
# options
'point': True/ False
'line': True/ False
'surf': True/ False
}
# beam cross-section view optinos
'cs': True,
},
}
fea.PreviewAnalysis(previewset)
Solver Set¶
# ==============================================================================
# solver set
# ==============================================================================
solverset = {"STEPSET":
{
'type': str()
# options
'table', # tabulated values, e.g. [step0, step1, ..., stepN]
'mode', # vibration modes
'freq', # frequency scale, np.linspace
'time', # time scale, np.linspace
'start': int(),
'end': int(),
'step': int(),
},
'SYMM': bool(), # matrix symmetric assembler
# options
True
False
# [ adv ]
'IBZ': np.array( # Brillouin zone
[[0.0, 0.0], # O
[np.pi, 0.0], # A
[np.pi, np.pi],# B
[0.0, 0.0]]), # O
# (B)
# /^
# / |
# / |
# / |
# / |
# / |
# </ |
# (O)--->(A)
# [ dev ]
'MP': bool(), # multi processing
# options
True
False
}
solverdata = fea.Solve(solverset)
Post-Process¶
# ==============================================================================
# post-process set
# ==============================================================================
postprocset = {
"SOLVERDATA": solverdata,
"COMPUTER": {
'structural': {
'displ': True/ False,
'stress': True/ False}
'thermal': {
'temp': True/ False,
'heatflux': True/ False}
# [dev]
'fluid',
},
"PLOTSET": {
'show': True/ False,
'filename': str(),
'savepng': True/ False},
"OUTPUT": {
'log': True/ False,
'get':{
'nelem': True/ False, # number of elements in the mesh
'nnode': True/ False, # number of nodes in the mesh
'inci': True/ False, # list of elements
'coord':True/ False, # list of nodes coodinate
'tabmat':True/ False, # list of material property
'tabgeo':True/ False, # list of geometry property
'bc_list':True/ False, # list of constraints (boundary conditions)
'lo_list':True/ False, # list of loads
'u_list':True/ False, # list of solutions
'numpy_decimals': int(), # decimals to print, e.g. 8
}
}}
postprocdata = fea.PostProcess(postprocset)
Tutorials¶
The main objective of these tutorials is to help the user explore the various modeling and simulation options that the project allows.
This tutorial was developed during the initial implementation of the project; therefore, the information contained herein is not 100% up-to-date with the latest version of myfempy. See the API documentation for a more modern version.
Mesh¶
Manual Mesh¶
'''
Geração da malha manual
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
force = {
'DEF': 'forceedge',
'DOF': 'fx',
'DIR': 'edgex',
'LOC': {'x': 0, 'y': 999, 'z': 0},
'VAL': [1.0],
}
bondcond = {
'DEF': 'fixed',
'DOF': 'all',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 0, 'z': 0},
}
# Retangulo estado plano (malha quad 4) 100 x 50 mm
elementos = [[1,"plane41","material","geometria",[1, 2, 3, 4]]]
coordenadas = [[1, 0, 0, 0],
[2, 100, 0, 0],
[3, 100, 50, 0],
[4, 0, 50, 0]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat],
"PROPGEO": [geo],
# "FORCES": [force],
# "BOUNDCOND": [bondcond],
# "QUADRATURE":{'meth':'gaussian','npp':8},
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01a', 'show': True, 'scale': 1, 'savepng': True, 'lines': True},
}
preview_plot(previewset, modelinfo)
# Retangulo estado plano (malha tria 4) 50 x 50 mm/ teste de qualidade 'Aspect Ratio'
elementos = [[1,"plane31","material","geometria",[1, 2, 4]],
[2,"plane31","material","geometria",[2, 3, 4]]]
coordenadas = [[1, 0, 0, 0],
[2, 50, 0, 0],
[3, 50, 40, 0],
[4, 0, 50, 0]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat],
"PROPGEO": [geo],
# "FORCES": [force],
# "BOUNDCOND": [bondcond],
# "QUADRATURE":{'meth':'gaussian','npp':8},
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01b', 'show': True, 'scale': 1, 'savepng': True, 'lines': True},
'QUALITY': {'show': True, 'method': 1, 'scale': 0.1},
}
preview_plot(previewset, modelinfo)
# Solido 100 x 100 x 100 mm
mat2 = {
"NAME": "material2",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'solid'
}
geo2 = {
"NAME": "geometria2",
"THICKN": 0.0
}
elementos = [[1,"solid81","material2","geometria2",[1, 2, 3, 4, 5, 6, 7, 8]]]
coordenadas = [[1, 0, 0, 0],
[2, 100, 0, 0],
[3, 100, 100, 0],
[4, 0, 100, 0],
[5, 0, 0, 100],
[6, 100, 0, 100],
[7, 100, 100, 100],
[8, 0, 100, 100]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat2],
"PROPGEO": [geo2],
# "FORCES": [force],
# "BOUNDCOND": [bondcond],
# "QUADRATURE":{'meth':'gaussian','npp':8},
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01c', 'show': True, 'scale': 1, 'savepng': True, 'lines': True},
'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
Legacy Mesh¶
'''
Geração da malha por meio da opção legacy
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
from math import pi
import sys
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'axial'
}
geo = {
"NAME": "geo",
"SEC":"I",
'DIM':{'b':100,'h':150,'t':5,'d':5}}
# Linha (malha beam 2) 500 mm
meshdata = {"LEGACY": {'lx': 400, 'mesh': 'line2', 'elem': 'beam21', 'nx': 2},
"PROPMAT": [mat],
"PROPGEO": [geo],
# "FORCES": [force],
# "BOUNDCOND": [bondcond],
# "QUADRATURE": {'meth': 'gaussian', 'npp': 4},
# "DOMAIN":'structural'
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_02a', 'show': True, 'scale': 1, 'savepng': True, 'lines': True, 'cs': True},
}
preview_plot(previewset, modelinfo)
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# Retangulo estado plano (malha quad 4) 100 x 50 mm
meshdata = {"LEGACY": {'lx': 100, 'ly': 50, 'mesh': 'quad4', 'elem': 'plane41', 'nx': 10, 'ny': 5},
"PROPMAT": [mat],
"PROPGEO": [geo],
# "FORCES": [force],
# "BOUNDCOND": [bondcond],
# "QUADRATURE": {'meth': 'gaussian', 'npp': 4},
# "DOMAIN":'structural'
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_02b', 'show': True, 'scale': 1, 'savepng': True, 'lines': True},
}
preview_plot(previewset, modelinfo)
# Retangulo estado plano (malha tria 3) 100 x 50 mm
meshdata = {"LEGACY": {'lx': 100, 'ly': 50, 'mesh': 'tria3', 'elem': 'plane31', 'nx': 10, 'ny': 5},
"PROPMAT": [mat],
"PROPGEO": [geo],
# "FORCES": [force],
# "BOUNDCOND": [bondcond],
# "QUADRATURE": {'meth': 'gaussian', 'npp': 4},
# "DOMAIN":'structural'
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_02c', 'show': True, 'scale': 1, 'savepng': True, 'lines': True},
}
preview_plot(previewset, modelinfo)
GMSH Mesh Basic¶
'''
Geração da malha com gmsh
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# Retangulo estado plano (malha quad 4) 100 x 50 mm
points = [[0, 0, 0],
[100, 0, 0],
[100, 50, 0],
[0, 50, 0]]
lines = [[1, 2],
[2, 3],
[3, 4],
[4, 1]]
plane = [[1, 2, 3, 4]]
meshdata = {"GMSH": {'filename': 'tutorial_03',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'tria3', #quad4
'elem': 'plane31', #plane41
'sizeelement': 5,
'meshmap': {'on': True,
'edge': [1, 2], #'all'
"numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_03',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {'point': True,
'edge': True}},
}
preview_plot(previewset, modelinfo)
# Placa L estado plano (malha quad 4) 100 x 100 mm
points = [[0, 0, 0],
[100, 0, 0],
[100, 50, 0],
[50, 50, 0],
[50, 100, 0],
[0, 100, 0],
]
lines = [[1, 2],
[2, 3],
[3, 4],
[4, 5],
[5, 6],
[6, 1],
]
plane = [[1, 2, 3, 4, 5, 6]]
meshdata = {"GMSH": {'filename': 'tutorial_03',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'quad4', #tria3
'elem': 'plane41', #plane31
'sizeelement': 5,
'meshmap': {'on': True,
'edge': 'all',
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_03',
'show': True,
'scale': 10,
'savepng': True,
'lines': True,
'plottags': {'point': True,
'edge': True}},
}
preview_plot(previewset, modelinfo)
# Placa H estado plano (malha quad 4) 100 x 100 mm
points = [
[0, 0, 0],
[20, 0, 0],
[20, 30, 0],
[80, 30, 0],
[80, 0, 0],
[100, 0, 0],
[100, 100, 0],
[80, 100, 0],
[80, 70, 0],
[20, 70, 0],
[20, 100, 0],
[0, 100, 0],
]
lines = [[1, 2],
[2, 3],
[3, 4],
[4, 5],
[5, 6],
[6, 7],
[7, 8],
[8, 9],
[9, 10],
[10, 11],
[11, 12],
[12, 1],
]
plane = [[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]]
meshdata = {"GMSH": {'filename': 'tutorial_03',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'quad4', #tria3
'elem': 'plane41', #plane31
'sizeelement': 5,
'meshmap': {'on': True,
'edge': 'all',
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_03',
'show': True,
'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
'point': True,
# 'edge': True
}
},
}
preview_plot(previewset, modelinfo)
GMSH Mesh Advanced¶
'''
Geração da malha com gmsh/ criação de curvas e furos e áreas compostas
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
import sys
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# 1/4 de placa com furo central (dia. 100 mm) estado plano (malha quad 4) 200 x 200 mm
points = [
[100, 0, 0],
[100, 100, 0],
[0, 100, 0],
]
circle = [[50, [0, 0, 0], ['0', 'Pi/2']]]
lines = [[1, 2],
[2, 3],
[3, 5],
[4, 1],
]
plane = [[1, 2, 3, 4, 5]]
meshdata = {"GMSH": {'filename': 'tutorial_04a',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'circle': circle,
'meshconfig': {
'mesh': 'quad4', #quad4
'elem': 'plane41', #plane41
'sizeelement': 5,
'meshmap': {'on': True,
'edge': 'all', #'all'
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04a',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {'point': True,
'edge': True}},
}
preview_plot(previewset, modelinfo)
# Placa completa multiplos furos (dia. 10, 10, 30 mm) estado plano (malha quad 4) 200 x 200 mm
points = [
[0, 0, 0],
[200, 0, 0],
[200, 200, 0],
[0, 200, 0],
]
lines = [[1, 2], # line 1
[2, 3], # line 2
[3, 4], # line 3
[4, 1], # line 4
]
circle = [[30, [150, 100, 0], ['0', '2*Pi']], # line 5
[10, [30, 150, 0], ['0', '2*Pi']], # line 6
[10, [30, 50, 0], ['0', '2*Pi']], # line 7
]
plane = [[1, 2, 3, 4],
[-5],
[-6],
[-7],
[5],
[6],
[7]
]
meshdata = {"GMSH": {'filename': 'tutorial_04b',
'pointlist': points,
'linelist': lines,
'circle': circle,
'planelist': plane,
'meshconfig': {
'mesh': 'quad4', #quad4
'elem': 'plane41', #plane41
'sizeelement': 10,
# 'meshmap': {'on': True,
# 'edge': [5], #'all'
# "numbernodes": 60,
# }
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04b',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
# 'edge': True,
'surf': True
}},
'LABELS': {'show': True, 'lines': True, 'scale': 0.2},
}
preview_plot(previewset, modelinfo)
# Placa completa furo quadrado central (dia. 50 mm) estado plano (malha quad 4) 200 x 200 mm
points = [
[0, 0, 0],
[200, 0, 0],
[200, 100, 0],
[0, 100, 0],
[75, 25, 0],
[125, 25, 0],
[125, 75, 0],
[75, 75, 0],
]
lines = [[1, 2], # line 1
[2, 3], # line 2
[3, 4], # line 3
[4, 1], # line 4
[5, 6], # line 5
[6, 7], # line 6
[7, 8], # line 7
[8, 5], # line 8
]
plane = [[1, 2, 3, 4],
[-5, -6, -7, -8],
# [5, 6, 7, 8]
]
meshdata = {"GMSH": {'filename': 'tutorial_04c',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'quad4', #quad4
'elem': 'plane41', #plane41
'sizeelement': 2,
'meshmap': {'on': True,
'edge': 'all', #'all'
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04c',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
'point': True,
# 'edge': True,
# 'surf': True
}},
}
preview_plot(previewset, modelinfo)
# Duas Placa composta estado plano (malha quad 4) 200 x 200 mm
points = [
[0, 0, 0], # point 1
[100, 0, 0], # point 2
[100, 100, 0], # point 3
[0, 100, 0], # point 4
[200, 0, 0], # point 5
[200, 100, 0], # point 6
]
lines = [[1, 2], # line 1
[2, 3], # line 2
[3, 4], # line 3
[4, 1], # line 4
[2, 5], # line 5
[5, 6], # line 6
[6, 3], # line 7
]
plane = [[1, 2, 3, 4], # plane 1
[5, 6, 7, 2], # plane 2
]
meshdata = {"GMSH": {'filename': 'tutorial_04d',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'quad4', #quad4
'elem': 'plane41', #plane41
'sizeelement': 5,
'meshmap': {'on': True,
'edge': [2, 6], #'all'
"numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04d',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
'edge': True,
# 'surf': True
}},
}
preview_plot(previewset, modelinfo)
# Duas Placa composta com furo em 0,0,0 (dia. 50 mm) estado plano (malha quad 4) 200 x 200 mm
points = [
[100, 0, 0], # point 1
[100, 100, 0], # point 2
[0, 100, 0], # point 3
[200, 0, 0], # point 4
[200, 100, 0], # point 5
]
lines = [
[6, 1], # line 1
[1, 2], # line 2
[2, 3], # line 3
[3, 7], # line 4
[1, 4], # line 5
[4, 5], # line 6
[5, 2], # line 7
]
circle = [[50, [0, 0, 0], ['0', 'Pi/2']], # line 8
]
plane = [[1, 2, 3, 4, 8], # plane 1
[5, 6, 7, 2], # plane 2
]
meshdata = {"GMSH": {'filename': 'tutorial_04e',
'pointlist': points,
'linelist': lines,
'circle': circle,
'planelist': plane,
'meshconfig': {
'mesh': 'tria3', #quad4
'elem': 'plane31', #plane41
'sizeelement': 10,
'meshmap': {'on': True,
'edge': [8], #'all'
"numbernodes": 20,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04e',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
'edge': True,
# 'surf': True
}},
}
preview_plot(previewset, modelinfo)
# disco (dia. 200 mm) com furo central (dia. 100 mm) estado plano (malha quad 4)
points = [
# [100, 0, 0],
# [100, 100, 0],
# [0, 100, 0],
]
lines = [
# [1, 2],
# [2, 3],
# [3, 5],
# [4, 1],
]
circle = [[200, [0, 0, 0], ['0', '2*Pi']],[100, [0, 0, 0], ['0', '2*Pi']]]
plane = [[1],[-2]]
meshdata = {"GMSH": {'filename': 'tutorial_04f',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'circle': circle,
'meshconfig': {
'mesh': 'tria3', #quad4
'elem': 'plane31', #plane41
'sizeelement': 40,
'meshmap': {'on': True,
'edge': 'all', #'all'
# "numbernodes": 80,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04f',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
'edge': True,
# 'surf': True
}},
# 'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
# 1/4 de placa com bara no furo central (dia. 100 mm) estado plano (malha quad 4) 200 x 200 mm
points = [
[100, 0, 0],
[100, 100, 0],
[0, 100, 0],
[0, 0, 0],
]
lines = [[1, 2],
[2, 3],
[3, 6],
[5, 1],
[4, 5],
[6, 4],
]
circle = [[50, [0, 0, 0], ['0', 'Pi/2']]]
plane = [[1, 2, 3, 4, 7],[5, 6, 7]]
meshdata = {"GMSH": {'filename': 'tutorial_04g',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'circle': circle,
'meshconfig': {
'mesh': 'quad4', #quad4
'elem': 'plane41', #plane41
'sizeelement': 10,
'meshmap': {'on': True,
'edge': [7], #'all'
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04g',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
# 'edge': True,
'surf': True,
}},
}
preview_plot(previewset, modelinfo)
#-----------------------------------------------------
# 1/4 de placa com furo central (dia. 100 mm) estado plano (malha quad 4) 200 x 200 mm malha mapeada
points = [
[100, 0, 0],
[100, 100, 0],
[0, 100, 0],
]
lines = [[1, 2], # 1
[2, 3], # 2
[3, 7], # 3
[4, 1], # 4
[2, 5], # 5
]
circle = [[50, [0, 0, 0], ['0', 'Pi/4']], # 6
[50, [0, 0, 0], ['Pi/4', 'Pi/2']]] # 7
plane = [[1, 5, 6, 4],
[5, 2, 3, 7]]
meshdata = {"GMSH": {'filename': 'tutorial_04h',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'circle': circle,
'meshconfig': {
'mesh': 'quad4', #quad4
'elem': 'plane41', #plane41
'sizeelement': 10,
'meshmap': {'on': True,
'edge': 'all', #'all'
"numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_04h',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
'point': True,
# 'edge': True
}},
# 'LABELS': {'show': True, 'lines': True, 'scale': 0.5},
}
preview_plot(previewset, modelinfo)
GMSH Mesh Solids¶
'''
Geração da malha com gmsh/ solidos
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
import sys
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# 1/4 de placa com furo central (dia. 100 mm) estado plano (malha quad 4) 200 x 200 mm, espessura da placa de 20 mm
points = [
[100, 0, 0],
[100, 100, 0],
[0, 100, 0],
]
lines = [[1, 2],
[2, 3],
[3, 5],
[4, 1],
]
arcs = [[50, [0, 0, 0], ['0', 'Pi/2']]]
plane = [[1, 2, 3, 4, 5]]
meshdata = {"GMSH": {'filename': 'tutorial_05a',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'arc': arcs,
'meshconfig': {
'mesh': 'tetr4', #tetr4 hexa8
'elem': 'solid41', #solid41 solid81
'sizeelement': 5,
'extrude': 20,
'meshmap': {'on': True,
'edge': 'all', #'all'
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_05a',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {'point': True,
'edge': True}},
}
preview_plot(previewset, modelinfo)
# Placa completa multiplos furos (dia. 10, 10, 30 mm) estado plano (malha quad 4) 200 x 200 mm, espessura da placa de 20 mm
points = [
[0, 0, 0],
[200, 0, 0],
[200, 200, 0],
[0, 200, 0],
]
lines = [[1, 2], # line 1
[2, 3], # line 2
[3, 4], # line 3
[4, 1], # line 4
]
arcs = [[30, [150, 100, 0], ['0', '2*Pi']], # line 5
[10, [30, 150, 0], ['0', '2*Pi']], # line 6
[10, [30, 50, 0], ['0', '2*Pi']], # line 7
]
plane = [[1, 2, 3, 4],
[-5],
[-6],
[-7],
# [5],
# [6],
# [7]
]
meshdata = {"GMSH": {'filename': 'tutorial_05b',
'pointlist': points,
'linelist': lines,
'arc': arcs,
'planelist': plane,
'meshconfig': {
'mesh': 'tetr4', #tetr4 hexa8
'elem': 'solid41', #solid41 solid81
'sizeelement': 5,
'extrude': 20,
'meshmap': {'on': True,
'edge': 'all', #'all'
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_05b',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
# 'edge': True,
'surf': True
}},
}
preview_plot(previewset, modelinfo)
# Placa completa furo quadrado central (dia. 50 mm) estado plano (malha quad 4) 200 x 200 mm, espessura da placa de 20 mm
points = [
[0, 0, 0],
[200, 0, 0],
[200, 100, 0],
[0, 100, 0],
[75, 25, 0],
[125, 25, 0],
[125, 75, 0],
[75, 75, 0],
]
lines = [[1, 2], # line 1
[2, 3], # line 2
[3, 4], # line 3
[4, 1], # line 4
[5, 6], # line 5
[6, 7], # line 6
[7, 8], # line 7
[8, 5], # line 8
]
plane = [[1, 2, 3, 4],
[-5, -6, -7, -8],
# [5, 6, 7, 8]
]
meshdata = {"GMSH": {'filename': 'tutorial_05c',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'hexa8', #tetr4 hexa8
'elem': 'solid81', #solid41 solid81
'sizeelement': 5,
'extrude': 20,
'meshmap': {'on': True,
'edge': 'all', #'all'
# "numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_05c',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
'point': True,
# 'edge': True,
# 'surf': True
}},
}
preview_plot(previewset, modelinfo)
# Duas Placa composta estado plano (malha quad 4) 200 x 200 mm, espessura da placa de 20 mm
points = [
[0, 0, 0], # point 1
[100, 0, 0], # point 2
[100, 100, 0], # point 3
[0, 100, 0], # point 4
[200, 0, 0], # point 5
[200, 100, 0], # point 6
]
lines = [[1, 2], # line 1
[2, 3], # line 2
[3, 4], # line 3
[4, 1], # line 4
[2, 5], # line 5
[5, 6], # line 6
[6, 3], # line 7
]
plane = [[1, 2, 3, 4], # plane 1
[5, 6, 7, 2], # plane 2
]
meshdata = {"GMSH": {'filename': 'tutorial_05d',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'tetr4', #tetr4 hexa8
'elem': 'solid41', #solid41 solid81
'sizeelement': 5,
'extrude': 20,
'meshmap': {'on': True,
'edge': [2, 6], #'all'
"numbernodes": 10,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_05d',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
'edge': True,
# 'surf': True
}},
}
preview_plot(previewset, modelinfo)
# ------------------------------------------------------
# Duas Placa composta com furo em 0,0,0 (dia. 50 mm) estado plano (malha quad 4) 200 x 200 mm, espessura da placa de 20 mm
points = [
[100, 0, 0], # point 1
[100, 100, 0], # point 2
[0, 100, 0], # point 3
[200, 0, 0], # point 4
[200, 100, 0], # point 5
]
lines = [
[6, 1], # line 1
[1, 2], # line 2
[2, 3], # line 3
[3, 7], # line 4
[1, 4], # line 5
[4, 5], # line 6
[5, 2], # line 7
]
arcs = [[50, [0, 0, 0], ['0', 'Pi/2']], # line 8
]
plane = [[1, 2, 3, 4, 8], # plane 1
[5, 6, 7, 2], # plane 2
]
meshdata = {"GMSH": {'filename': 'tutorial_05e',
'pointlist': points,
'linelist': lines,
'arc': arcs,
'planelist': plane,
'meshconfig': {
'mesh': 'tetr4', #tetr4 hexa8
'elem': 'solid41', #solid41 solid81
'sizeelement': 5,
'extrude': 20,
'meshmap': {'on': False,
'edge': 'all',
"numbernodes": 20,
}
}
},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_05e',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {
# 'point': True,
# 'edge': True,
'surf': True
}},
}
preview_plot(previewset, modelinfo)
GMSH Mesh Importing CAD models (.stp/.step)¶
'''
Geração da malha com gmsh/ importação de modelos cad (.stp/ .step)
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
import sys
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# import cad .stp
meshdata = {"GMSH": {'filename': 'tutorial_06a',
'cadimport': {'object': 'tutorial_06_cad.stp'},
'meshconfig': {'mesh': 'tetr4', 'elem': 'solid41', 'sizeelement': 10}},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_06a',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {'point': False,
'edge': False}},
}
preview_plot(previewset, modelinfo)
GMSH Mesh Importing mesh files (.msh2)¶
'''
Geração da malha com gmsh/ importação de malha (.vtk) externa
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
import sys
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# import mesh .msh2
meshdata = {"GMSH": {'filename': 'tutorial_06b',
'meshimport':{'object':'tutorial_06_mesh'},
'meshconfig': {'mesh': 'tetr4', 'elem': 'solid41', 'sizeelement': 10}},
"PROPMAT": [mat],
"PROPGEO": [geo],
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_06b',
'show': True, 'scale': 10,
'savepng': True,
'lines': True,
'plottags': {'point': False,
'surf': True}},
}
preview_plot(previewset, modelinfo)
Simulations¶
Structural Static¶
'''
Simulação estática com malha manual
'''
# Imports
from myfempy import ModelGen
from myfempy import Solver
from myfempy import PostProcess
from myfempy import postproc_plot
from myfempy import preview_plot
# Definição do material, geometria
mat = {
"NAME": "material",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo = {
"NAME": "geometria",
"THICKN": 1.0
}
# Retangulo estado plano (malha manual quad 4) 100 x 50 mm com forcas e bc. nodais
f1 = {
'DEF': 'forcenode',
'DOF': 'fy',
'DIR': 'node',
'LOC': {'x': 100, 'y': 0, 'z': 0},
'VAL': [-1000.0],
}
f2 = {
'DEF': 'forcenode',
'DOF': 'fy',
'DIR': 'node',
'LOC': {'x': 100, 'y': 50, 'z': 0},
'VAL': [-1000.0],
}
bc1 = {
'DEF': 'fixed',
'DOF': 'all',
'DIR': 'node',
'LOC': {'x': 0, 'y': 0, 'z': 0},
}
bc2 = {
'DEF': 'fixed',
'DOF': 'ux',
'DIR': 'node',
'LOC': {'x': 0, 'y': 50, 'z': 0},
}
elementos = [[1,"plane41","material","geometria",[1, 2, 3, 4]]]
coordenadas = [[1, 0, 0, 0],
[2, 100, 0, 0],
[3, 100, 50, 0],
[4, 0, 50, 0]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat],
"PROPGEO": [geo],
"FORCES": [f1,f2],
"BOUNDCOND": [bc1,bc2],
"QUADRATURE":{'meth':'gaussian','npp':4},
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01a', 'show': True, 'scale': 10, 'savepng': True, 'lines': True},
'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
# Retangulo estado plano (malha manual quad 4) 100 x 50 mm com forcas distribuida e bc. nodais
f1 = {
'DEF': 'forceedge',
'DOF': 'fy',
'DIR': 'edgex',
'LOC': {'x': 100, 'y': 999, 'z': 0},
'VAL': [-1000.0],
}
bc1 = {
'DEF': 'fixed',
'DOF': 'all',
'DIR': 'node',
'LOC': {'x': 0, 'y': 0, 'z': 0},
}
bc2 = {
'DEF': 'fixed',
'DOF': 'ux',
'DIR': 'node',
'LOC': {'x': 0, 'y': 50, 'z': 0},
}
elementos = [[1,"plane41","material","geometria",[1, 2, 3, 4]]]
coordenadas = [[1, 0, 0, 0],
[2, 100, 0, 0],
[3, 100, 50, 0],
[4, 0, 50, 0]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat],
"PROPGEO": [geo],
"FORCES": [f1],
"BOUNDCOND": [bc1,bc2],
"QUADRATURE":{'meth':'gaussian','npp':4},
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01b', 'show': True, 'scale': 10, 'savepng': True, 'lines': True},
'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
# Retangulo estado plano (malha legacy quad 4 10x5) 100 x 50 mm com forcas e bc distribuida
f1 = {
'DEF': 'forceedge',
'DOF': 'fy',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 50, 'z': 0},
'VAL': [-1000.0],
}
bc1 = {
'DEF': 'fixed',
'DOF': 'ux',
'DIR': 'edgex',
'LOC': {'x': 0, 'y': 999, 'z': 0},
}
bc2 = {
'DEF': 'fixed',
'DOF': 'all',
'DIR': 'node',
'LOC': {'x': 100, 'y': 0, 'z': 0},
}
meshdata = {"LEGACY": {'lx': 100, 'ly': 50, 'mesh': 'quad4', 'elem': 'plane41', 'nx': 10, 'ny': 5},
"PROPMAT": [mat],
"PROPGEO": [geo],
"FORCES": [f1],
"BOUNDCOND": [bc1,bc2],
"QUADRATURE": {'meth': 'gaussian', 'npp': 4},
# "DOMAIN":'structural'
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01c', 'show': True, 'scale': 10, 'savepng': True, 'lines': True},
'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
# Retangulo estado plano (malha manual quad 4) 100 x 50 mm com massa e mola concentrada
f1 = {
'DEF': 'forcenode',
'DOF': 'spring2ground',
'DIR': 'node',
'LOC': {'x': 100, 'y': 0, 'z': 0},
'VAL': [100.0],
}
f2 = {
'DEF': 'forcenode',
'DOF': 'damper2ground',
'DIR': 'node',
'LOC': {'x': 100, 'y': 50, 'z': 0},
'VAL': [100.0],
}
f3 = {
'DEF': 'forcenode',
'DOF': 'masspoint',
'DIR': 'node',
'LOC': {'x': 100, 'y': 50, 'z': 0},
'VAL': [1.0],
}
bc1 = {
'DEF': 'fixed',
'DOF': 'all',
'DIR': 'node',
'LOC': {'x': 0, 'y': 0, 'z': 0},
}
bc2 = {
'DEF': 'fixed',
'DOF': 'ux',
'DIR': 'node',
'LOC': {'x': 0, 'y': 50, 'z': 0},
}
elementos = [[1,"plane41","material","geometria",[1, 2, 3, 4]]]
coordenadas = [[1, 0, 0, 0],
[2, 100, 0, 0],
[3, 100, 50, 0],
[4, 0, 50, 0]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat],
"PROPGEO": [geo],
"FORCES": [f1,f2,f3],
"BOUNDCOND": [bc1,bc2],
"QUADRATURE":{'meth':'gaussian','npp':4},
}
modelinfo = ModelGen.get_model(meshdata)
previewset = {'RENDER': {'filename': 'tutorial_01a', 'show': True, 'scale': 10, 'savepng': True, 'lines': True},
# 'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
# Retangulo estado plano dois materiais (malha manual quad 4) 200 x 100 mm com forcas e bc. nodais
mat2 = {
"NAME": "material2",
"VXX": 0.25,
"EXX": 200E6,
"MAT": 'isotropic',
"DEF": 'planestress'
}
geo2 = {
"NAME": "geometria2",
"THICKN": 1.0
}
f1 = {
'DEF': 'forcenode',
'DOF': 'fy',
'DIR': 'node',
'LOC': {'x': 0, 'y': 100, 'z': 0},
'VAL': [-100.0],
}
f2 = {
'DEF': 'forcenode',
'DOF': 'fy',
'DIR': 'node',
'LOC': {'x': 100, 'y': 100, 'z': 0},
'VAL': [-100.0],
}
f3 = {
'DEF': 'forcenode',
'DOF': 'fy',
'DIR': 'node',
'LOC': {'x': 200, 'y': 100, 'z': 0},
'VAL': [-100.0],
}
bc1 = {
'DEF': 'fixed',
'DOF': 'all',
'DIR': 'node',
'LOC': {'x': 200, 'y': 0, 'z': 0},
}
bc2 = {
'DEF': 'fixed',
'DOF': 'ux',
'DIR': 'node',
'LOC': {'x': 0, 'y': 0, 'z': 0},
}
bc3 = {
'DEF': 'fixed',
'DOF': 'ux',
'DIR': 'node',
'LOC': {'x': 0, 'y': 100, 'z': 0},
}
elementos = [[1,"plane41","material","geometria",[1, 2, 3, 4]],
[2,"plane41","material2","geometria2",[2, 5, 6, 3]]]
coordenadas = [[1, 0, 0, 0],
[2, 100, 0, 0],
[3, 100, 100, 0],
[4, 0, 100, 0],
[5, 200, 0, 0],
[6, 200, 100, 0]]
meshdata = {"ADD121": elementos,
"NODELIST": coordenadas,
"PROPMAT": [mat,mat2],
"PROPGEO": [geo,geo2],
"FORCES": [f1,f2,f3],
"BOUNDCOND": [bc1,bc2,bc3],
"QUADRATURE":{'meth':'gaussian','npp':4},
}
modelinfo = ModelGen.get_model(meshdata)
print(modelinfo['inci'])
previewset = {'RENDER': {'filename': 'tutorial_01a', 'show': True, 'scale': 10, 'savepng': True, 'lines': True},
'LABELS': {'show': True, 'lines': True, 'scale': 10},
}
preview_plot(previewset, modelinfo)
Structural Static Solid Gmsh API¶
'''
Simulação estática em sólido com geração de malha com gmsh
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
# from os import environ
# environ['OMP_NUM_THREADS'] = '1' #win
import numpy as np
import matplotlib.pyplot as plt
from myfempy import newAnalysis
from myfempy import SteadyStateLinear, SteadyStateLinearIterative
from time import time
# ===============================================================================
# FEA
# ===============================================================================
fea = newAnalysis(SteadyStateLinear) # StaticLinearIterative FASTER THEN StaticLinear
# MODEL SET
mat1 = {
"NAME": "aco",
"VXY": 0.33,
"EXX": 71E3,
# "RHO": 1, # kg/mm^2
}
geo2 = {"NAME": "secao1",}
# MODEL SET
# mesh tetra 4
nodes = [[1, 0, 0, 0],
[2, 1, 0, 0],
[3, 0, 1, 0],
[4, 0, 0, 1]]
# nodes = [[1, 500., 500., 0.],
# [2, 0., 500., 250.],
# [3, 0., 1000., 0.],
# [4, 0., 0., 0.]]
conec = [[1, 1, 1, 1, 2, 3, 4]]
# mesh hexa 8
nodes = [
[1, 0, 0, 0], # Nó 1
[2, 1, 0, 0], # Nó 2
[3, 1, 1, 0], # Nó 3
[4, 0, 1, 0], # Nó 4
[5, 0, 0, 1], # Nó 5
[6, 1, 0, 1], # Nó 6
[7, 1, 1, 1], # Nó 7
[8, 0, 1, 1] # Nó 8
]
conec = [[1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 8]]
# # gmsh config
LX = 200
LY = 20
LZ = 20
points = [
[0, 0, 0],
[LX, 0, 0],
[LX, LY, 0],
[0, LY, 0],
]
lines = [[1, 2],
[2, 3],
[3, 4],
[4, 1],
]
plane = [[1, 2, 3, 4]]
modeldata = {
# "MESH": {
# 'TYPE': 'manual',
# 'COORD': nodes,
# 'INCI': conec,
# },
"MESH": {
'TYPE': 'gmsh',
'filename': 'data_mesh',
'pointlist': points,
'linelist': lines,
'planelist': plane,
'meshconfig': {
'mesh': 'hexa8',
# "numbernodes": 21,
'sizeelement': 10,
'extrude': LZ,
'meshmap': {'on': True,
'edge': 'all',
"numbernodes": [20],
}
}
},
"ELEMENT": {
'TYPE': 'structsolid',
'SHAPE': 'hexa8',
# 'INTGAUSS': 4,
},
"MATERIAL": {
"MAT": 'solidelastic',
"TYPE": 'isotropic',
"PROPMAT": [mat1],
},
"GEOMETRY": {
"GEO": 'solid',
"PROPGEO": [geo2],
},
}
fea.Model(modeldata)
inci = fea.getInci()
coord = fea.getCoord()
# tabmat = fea.getTabmat()
# tabgeo = fea.getTabgeo()
# intgauss = fea.getIntGauss()
# print(intgauss)
# vol = fea.getElementVolume(inci, coord, tabgeo)
# print('Volume ',np.sum(vol))
# reg = fea.getRegions()
# print(reg)
# sys.exit()
# element_number = 0
# st = time()
# ke = fea.getElemStifLinearMat(inci, coord, tabmat, tabgeo, intgauss, element_number)
# print('time KE ', time() - st)
# print(np.array2string(ke, separator=', '))
# print(np.allclose(ke, ke.T, rtol=1e-05, atol=1e-08))
# print(np.all(np.linalg.eigvals(ke)))
# st = time()
# kg = fea.getGlobalMatrix(inci, coord, tabmat, tabgeo, intgauss, SYMM=False)
# t_kg = time()-st
# print('time numpy sym ', t_kg)
# kg = kg['stiffness'].todense()
# print(kg.todense())
# print('kg is sym ', scipy.linalg.issymmetric(kg, atol = 1e-09))
# plt.figure(2)
# plt.spy(kg, markersize=4)
# plt.show()
# sys.exit()
f1 = {
'TYPE': 'forcenode',
'DOF': 'fx',
'DIR': 'node',
'LOC': {'x': 0, 'y': 0, 'z': 0},
'VAL': [1],
}
f1 = {
'TYPE': 'forcenode',
'DOF': 'fx',
'DIR': 'surfyz',
'LOC': {'x': LX, 'y': 999, 'z': 999},
'VAL': [1],
}
f2 = {
'TYPE': 'forcesurf',
'DOF': 'fx',
'DIR': 'surfyz',
'LOC': {'x': LX, 'y': 999, 'z': 999},
'VAL': [1],
}
f2b = {
'TYPE': 'forcesurf',
'DOF': 'fy',
'DIR': 'plane',
'TAG': 4,
'VAL': [-1.0],
}
f3 = {
'TYPE': 'forcesurf',
'DOF': 'pressure',
'DIR': 'surfyz',
'LOC': {'x': LX, 'y': 999, 'z': 999},
'VAL': [1],
}
f3pt = {
'TYPE': 'forcesurf',
'DOF': 'pressure',
'DIR': 'plane',
'TAG': 4,
'VAL': [-1.0],
}
bc_node = {
'TYPE': 'fixed',
'DOF': 'full',
'DIR': 'node',
'LOC': {'x': 1, 'y': 1, 'z': 1},
}
bc = {
'TYPE': 'fixed',
'DOF': 'full',
'DIR': 'surfyz',
'LOC': {'x': 0, 'y': 999, 'z': 999},
}
physicdata = {
"PHYSIC": {"DOMAIN": "structural", # 'fluid' 'thermal'; "COUPLING": 'fsi'
"LOAD": [f2b],
"BOUNDCOND": [bc],
},
}
fea.Physic(physicdata)
# regions = fea.getRegions()
# print(regions)
# loadaply = fea.getLoadApply()
# # bcaply = fea.getBCApply()
# print(loadaply)
# # print(bcaply)
# # print(fea.getLoadArray(loadaply))
# # print(fea.getDirichletNH(bcaply))
# print('forca total [N]', np.sum(loadaply[:,2]))
# sys.exit()
previewset = {'RENDER': {'filename': 'beam', 'show': True, 'scale': 2, 'savepng': True, 'lines': True,
'plottags': {'point': True}
},
# 'LABELS': {'show': True, 'lines': True, 'scale': 1},
}
fea.PreviewAnalysis(previewset)
# sys.exit()
# #-------------------------------- SOLVER -------------------------------------#
solverset = {"STEPSET": {'type': 'table', # mode, freq, time ...
'start': 0,
'end': 1,
'step': 1},
'SYMM':False,
# 'MP':True,
}
solverdata = fea.Solve(solverset)
print(max(abs(solverdata['solution']['U'])))
postprocset = {"SOLVERDATA": solverdata,
"COMPUTER": {'structural': {'displ': True, 'stress': True}},
"PLOTSET": {'show': True, 'filename': 'test_shakedown', 'savepng': True},
# "TRACKER": {'point': {'x': 0, 'y': 0, 'z': 0, 'dof':1}},
"OUTPUT": {'log': True, 'get': {'nelem': True, 'nnode': True}},
}
postprocdata = fea.PostProcess(postprocset)
# print(postprocdata["solution"])
Structural Vibration¶
'''
Simulação de vibração modal
'''
import sys
# setting path
sys.path.append('../myfempy')
import sys
import numpy as np
import matplotlib.pyplot as plt
from myfempy import newAnalysis
from myfempy import DynamicEigenLinear
import time
#===============================================================================
# FEA
#===============================================================================
fea = newAnalysis(DynamicEigenLinear)
mat1 = {
"NAME": "mat1",
"VXY": 0.3,
"EXX": 210E6,
"RHO": 7850,
}
mat2 = {
"NAME": "mat2",
"VXX": 0.3,
"EXX": 1E6,
"RHO": 1,
}
geo = {
"NAME": "espessura1",
"THICKN": 0.1,
# "DIM": [b, h, t, d],
}
# MODEL SET
# nelx = 10
# nely = 5
# nodes = [[1, 0, 0, 0],
# [2, nelx, 0, 0],
# [3, nelx, nely, 0],
# [4, 0, nely, 0],
# # [5, 1, 1, 0],
# # [6, 0, 1, 0],
# ]
# conec = [[1, 1, 1, 1, 2, 3, 4],
# # [2, 2, 1, 2, 3, 4, 5],
# ]
# MODEL SET
LX = 8 # 80
LY = 8 # 60
nelx = 4 # 40 # 80 # 160
nely = 4 # 30 # 60 # 120
# esize = 0.25
# # gmsh config
# points = [
# [0, 0, 0],
# [LX, 0, 0],
# [LX, LY, 0],
# [0, LY, 0]
# ]
# lines = [[1, 2],
# [2, 3],
# [3, 4],
# [4, 1],
# ]
# plane = [[1, 2, 3, 4]]
modeldata = {
# "MESH": {
# 'TYPE': 'add',
# 'COORD': nodes,
# 'INCI': conec,
# },
"MESH": {
'TYPE': 'legacy',
'LX': LX,
'LY': LY,
'NX': nelx,
'NY': nely,
},
# "MESH": {
# 'TYPE': 'gmsh',
# 'filename': 'test_mfoop_gmsh_01',
# # "meshimport": 'object_dir',
# 'pointlist': points,
# 'linelist': lines,
# 'planelist': plane,
# # 'arc': arcs,
# 'meshconfig': {
# 'mesh': 'quad4', #quad4
# 'sizeelement': 2 * esize,
# # 'extrude': 20,
# 'meshmap': {'on': True,
# 'edge': 'all', #'all'
# # "numbernodes": 10,
# }}
# },
"ELEMENT": {
'TYPE': 'structplane',
'SHAPE': 'quad4',
'INTGAUSS': 4,
},
"MATERIAL": {
"MAT": 'planestress',
"TYPE": 'isotropic',
"PROPMAT": [mat2],
},
"GEOMETRY": {
"GEO": 'thickness',
"PROPGEO": [geo],
},
}
fea.Model(modeldata)
forcespring_left = {
'TYPE': 'forcenode',
'DOF': 'spring2ground',
'DIR': 'point',
'TAG': 1,
'VAL': [1000.0],
}
forcespring_right = {
'TYPE': 'forcenode',
'DOF': 'spring2ground',
'DIR': 'point',
'TAG': 2,
'VAL': [1000.0],
}
massadd = {
'TYPE': 'forcenode',
'DOF': 'masspoint',
'DIR': 'node', #'point', #'node',
'LOC': {'x': LX/2, 'y': LY/2, 'z': 0},
'VAL': [1.4],
}
bc_node_left = {
'TYPE': 'fixed',
'DOF': 'full',
'DIR': 'node', # node
'LOC': {'x': 0, 'y': LY/2, 'z': 0},
}
bc_node_right = {
'TYPE': 'fixed',
'DOF': 'full',
'DIR': 'node',
'LOC': {'x': LX, 'y': LY/2, 'z': 0},
}
physicdata = {
"PHYSIC": {"DOMAIN": "structural",
"LOAD": [],
"BOUNDCOND": [],
},
}
fea.Physic(physicdata)
# loadaply = fea.getLoadApply()
# # bcaply = fea.getBCApply()
# print(loadaply)
# print(bcaply)
# print(fea.getLoadArray(loadaply))
# print(fea.getDirichletNH(bcaply))
# print('forca total [N]', np.sum(loadaply[:,2]))
previewset = {'RENDER': {'filename': 'squad', 'show': True, 'scale': 10, 'savepng': True, 'lines': False,
# 'plottags': {
# # 'line': True
# 'point': True
# }
},
# 'LABELS': {'show': True, 'lines': True, 'scale': 1},
}
fea.PreviewAnalysis(previewset)
# sys.exit()
# #-------------------------------- SOLVER -------------------------------------#
solverset = {"STEPSET": {'type': 'mode', # mode, freq, time ...
'start': 0,
'end': 10,
'step': 1},
'SYMM':False,
# 'MP':True,
}
solverdata = fea.Solve(solverset)
print(solverdata['solution']['FREQ'])
postprocdata = {"SOLVERDATA": solverdata,
"COMPUTER": {
'structural': {'modes': True},
# 'structural': {'frf': True}
},
"PLOTSET": {'filename': 'test_dynamic', 'savepng': True},
# "TRACKER": {'frf': {'x': 16, 'y': 1, 'z': 0, 'dof':0}},
"OUTPUT": {'log': True, 'get': {'nelem': True, 'nnode': True}},
}
postporc_result = fea.PostProcess(postprocdata)
Heat¶
'''
Simulação com tranferencia de calor em sólidos
Necessario instalação prévia do gmsh (não nativo do myfempy)
'''
from myfempy import newAnalysis
from myfempy import SteadyStateLinear, SteadyStateLinearIterative
from time import time
# ===============================================================================
# FEA
# ===============================================================================
fea = newAnalysis(SteadyStateLinearIterative)
mat1 = {
"NAME": "aco",
"KXX": 0.0605, # W/mm·°C
"KYY": 0.0605,
}
geo = {
"NAME": "espessura",
"THICKN": 0.5, # mm
# "DIM": [b, h, t, d],
}
# MODEL SET
LX = 20 # 80 # mm
LY = 20 # 60
nelx = 80 # 40 # 80 # 160
nely = 80 # 30 # 60 # 120
# esize = 1 #LX/nelx
# MODEL SET
# example 10.4 Logan
# nodes = [[1, 0, 0, 0],
# [2, LX, 0, 0],
# [3, LX, LY, 0],
# [4, 0, LY, 0]]
# conec = [[1, 1, 1, 1, 2, 3, 4]]
# gmsh config
points = [
[0, 0, 0],
[20, 0, 0],
[20, 20, 0],
[0, 20, 0]
]
lines = [[1, 2],
[2, 3],
[3, 4],
[4, 1],
]
plane = [[1, 2, 3, 4]]
modeldata = {
# "MESH": {
# 'TYPE': 'add',
# 'COORD': nodes,
# 'INCI': conec,
# },
# "MESH": {
# 'TYPE': 'legacy',
# 'LX': LX,
# 'LY': LY,
# 'NX': nelx,
# 'NY': nely,
# },
"MESH": {
'TYPE': 'gmsh',
'filename': 'test_line_force_x',
# "meshimport": 'object_dir',
'pointlist': points,
'linelist': lines,
'planelist': plane,
# 'arc': arcs,
'meshconfig': {
'mesh': 'quad4', #quad4
'sizeelement': 2,
# 'extrude': 5,
'meshmap': {'on': True,
'edge': 'all', #'all'
"numbernodes": 4,
}
}
},
"ELEMENT": {
'TYPE': 'heatplane',
'SHAPE': 'quad4',
# 'INTGAUSS': 8,
},
"MATERIAL": {
"MAT": 'heatplane',
"TYPE": 'isotropic',
"PROPMAT": [mat1],
},
"GEOMETRY": {
"GEO": 'thickness',
"PROPGEO": [geo],
},
}
fea.Model(modeldata)
# inci = fea.getInci()
# coord = fea.getCoord()
# tabmat = fea.getTabmat()
# tabgeo = fea.getTabgeo()
# intgauss = fea.getIntGauss()
# element_number = 0
# st = time()
# ke = fea.getElemStifLinearMat(inci, coord, tabmat, tabgeo, intgauss, element_number)
# print('time KE ', time() - st)
# print(np.array2string(ke, separator=', '))
# # sys.exit()
# kg = fea.getGlobalMatrix(inci, coord, tabmat, tabgeo, intgauss)
# kg = kg['stiffness']
# print(kg.todense())
# import matplotlib.pylab as plt
# import scipy
# print('kg is sym ', scipy.linalg.issymmetric(kg.todense(), atol = 1e-09))
# # print('kg_cy is close to kg ',np.allclose(kg.todense(), kg_cy.todense()))
# plt.figure(2)
# plt.spy(kg.todense(), markersize=4)
# plt.show()
# sys.exit()
q = 1.0 #1
Q = 0.2 # W/mm3
h = 0.000005 # W/mm2.ºC
Tinf = 1000
heatflux = {
'TYPE': 'heatfluxedge',
'DOF': 'heatflux',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 0, 'z': 0},
'VAL': [q],
}
heatgen = {
'TYPE': 'heatgeneration',
'DOF': 'heatflux',
'DIR': 'node',
# 'LOC': {'x': 0, 'y': 999, 'z': 0},
'VAL': [Q],
}
convc_s1 = {
'TYPE': 'convectionedge',
'DOF': 'convection',
'DIR': 'edgex',
'LOC': {'x': LX, 'y': 999, 'z': 0},
'VAL': [h*Tinf],
}
convc_s2 = {
'TYPE': 'convectionedge',
'DOF': 'convection',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 0, 'z': 0},
'VAL': [h*Tinf],
}
convc_s3 = {
'TYPE': 'convectionedge',
'DOF': 'convection',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': LY, 'z': 0},
'VAL': [h*Tinf],
}
bc1 = {
'TYPE': 'insulated',
'DOF': 'full',
'DIR': 'edgex',
'LOC': {'x': 0, 'y': 999, 'z': 0},
}
bc2 = {
'TYPE': 'insulated',
'DOF': 'full',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 0, 'z': 0},
}
bc3 = {
'TYPE': 'insulated',
'DOF': 'full',
'DIR': 'edgex',
'LOC': {'x': LX, 'y': 999, 'z': 0},
}
bcNH1 = {
'TYPE': 'temperature',
'DOF': 't',
'DIR': 'edgex',
'LOC': {'x': 0, 'y': 999, 'z': 0},
'VAL': [180.0],
}
bcNH2 = {
'TYPE': 'temperature',
'DOF': 't',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': LY, 'z': 0},
'VAL': [1.0],
}
physicdata = {
"PHYSIC": {"DOMAIN": "thermal", # 'fluid' 'thermal'; "COUPLING": 'fsi'
"LOAD": [convc_s1, heatflux],
"BOUNDCOND": [bcNH2],
}
}
fea.Physic(physicdata)
loadaply = fea.getLoadApply()
bcaply = fea.getBCApply()
print(loadaply)
print(bcaply)
previewset = {'RENDER': {'filename': 'heat_transfer', 'show': True, 'scale': 10, 'savepng': True, 'lines': False,
# 'plottags': {'line': True}
},
# 'LABELS': {'show': True, 'lines': True, 'scale': 1},
}
fea.PreviewAnalysis(previewset)
# nodes = fea.getNodesFromRegions(1, 'plane')
# print(nodes)
# elem = fea.getElementFromNodesList(nodes)
# print(elem)
# #-------------------------------- SOLVER -------------------------------------#
solverset = {"STEPSET": {'type': 'table', # mode, freq, time ...
'start': 0,
'end': 1,
'step': 1},
'SYMM':True,
# 'MP':True,
}
solverdata = fea.Solve(solverset)
postprocset = {"SOLVERDATA": solverdata,
"COMPUTER": {'thermal': {'temp': True, 'heatflux': True}},
"PLOTSET": {'show': True, 'filename': 'test_shakedown', 'savepng': True},
# "TRACKER": {'point': {'x': 0, 'y': 0, 'z': 0, 'dof':1}},
"REPORT": {'log': True, 'get':{
'nelem': True,
'nnode': True,
'inci': True,
'coord':True,
'tabmat':True,
'tabgeo':True,
'bc_list':True,
'lo_list':True,
'u_list': True,
}
}}
postprocdata = fea.PostProcess(postprocset)
Thermo Mechanical Coupling¶
'''
Simulação com iteração termo-mecânico em sólido 3D
'''
from os import environ
environ['OMP_NUM_THREADS'] = '1' #win
import numpy as np
import matplotlib.pyplot as plt
from myfempy import newAnalysis
from myfempy import SteadyStateLinear, SteadyStateLinearIterative
from time import time
# ===============================================================================
# FEA
# ===============================================================================
fea = newAnalysis(SteadyStateLinear)
steel = {
"NAME": "steel",
"KXX": 0.0605, # W/mm·°C
"KYY": 0.0605,
"KZZ": 0.0605,
"VXY": 0.30,
"EXX": 200E3, # N/mm^2 --> MPa
"RHO": 7.85E-06, # kg/mm^3
"CTE": 12E-6 # (mm/mm) /°C
}
copper = {
"NAME": "copper",
"KXX": 0.4230, # W/mm·°C
"KYY": 0.4230,
"KZZ": 0.4230,
"VXY": 0.33,
"EXX": 110E3, # N/mm^2 --> MPa
"RHO": 7.85E-06, # kg/mm^3
"CTE": 17E-6 # (mm/mm) /°C
}
geo = {
"NAME": "espessura",
"THICKN": 5, # mm
}
points = [
[0, 0, 0],
[200, 0, 0],
[0, 10, 0],
[200, 10, 0],
[0, 20, 0],
[200, 20, 0]
]
lines = [[1, 2],
[2, 4],
[3, 4],
[3, 1],
[6, 4],
[5, 6],
[5, 3],
]
plane = [[1, 2, 3, 4], [3, 5, 6, 7]]
modeldata = {
"MESH": {
'TYPE': 'gmsh',
'filename': 'test_line_force_x',
# "meshimport": 'object_dir',
'pointlist': points,
'linelist': lines,
'planelist': plane,
# 'arc': arcs,
'meshconfig': {
'mesh': 'hexa8', #quad4
'sizeelement': 4,
'extrude': 20,
'meshmap': {'on': True,
'edge': 'all', #'all'
"numbernodes": [10],
}
}
},
"ELEMENT": {
'TYPE': 'heatsolid',
'SHAPE': 'hexa8',
# 'INTGAUSS': 7,
},
"MATERIAL": {
"MAT": 'heatsolid',
"TYPE": 'isotropic',
"PROPMAT": [copper, steel],
},
"GEOMETRY": {
"GEO": 'thickness',
"PROPGEO": [geo, geo],
},
}
fea.Model(modeldata)
# inci = fea.getInci()
# coord = fea.getCoord()
# tabmat = fea.getTabmat()
# tabgeo = fea.getTabgeo()
# intgauss = fea.getIntGauss()
# vol = fea.getElementVolume(inci, coord, tabgeo, intgauss)
# print(np.sum(vol))
# # sys.exit()
q = -0.1 #1 # W/mm2
# Q = 0.2 # W/mm3
h = 0.000005 # W/mm2.ºC 0.000005 (air)
Tinf = 20 # ºC
heatflux = {
'TYPE': 'heatfluxsurf',
'DOF': 'heatflux',
'DIR': 'edgex',
'LOC': {'x': 200, 'y': 999, 'z': 0},
'VAL': [0*q],
}
# heatgen = {
# 'TYPE': 'heatgeneration',
# 'DOF': 'heatflux',
# 'DIR': 'node',
# # 'LOC': {'x': 0, 'y': 999, 'z': 0},
# 'VAL': [Q],
# }
convc_s1 = {
'TYPE': 'convectionedge',
'DOF': 'convection',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 0, 'z': 0},
'VAL': [-h*Tinf],
}
convc_s2 = {
'TYPE': 'convectionedge',
'DOF': 'convection',
'DIR': 'edgey',
'LOC': {'x': 999, 'y': 20, 'z': 0},
'VAL': [-h*Tinf],
}
# convc_s3 = {
# 'TYPE': 'convectionedge',
# 'DOF': 'convection',
# 'DIR': 'edgey',
# 'LOC': {'x': 999, 'y': LY, 'z': 0},
# 'VAL': [h*Tinf],
# }
# bc1 = {
# 'TYPE': 'insulated',
# 'DOF': 'full',
# 'DIR': 'edgex',
# 'LOC': {'x': 0, 'y': 999, 'z': 0},
# }
# bc2 = {
# 'TYPE': 'insulated',
# 'DOF': 'full',
# 'DIR': 'edgey',
# 'LOC': {'x': 999, 'y': 0, 'z': 0},
# }
# bc3 = {
# 'TYPE': 'insulated',
# 'DOF': 'full',
# 'DIR': 'edgex',
# 'LOC': {'x': LX, 'y': 999, 'z': 0},
# }
bcNH1 = {
'TYPE': 'temperature',
'DOF': 't',
'DIR': 'surfyz',
'LOC': {'x': 0, 'y': 999, 'z': 999},
'VAL': [200.0],
}
bcNH2 = {
'TYPE': 'temperature',
'DOF': 't',
'DIR': 'edgex',
'LOC': {'x': 200, 'y': 999, 'z': 0},
'VAL': [200.0],
}
# bcNH3 = {
# 'TYPE': 'temperature',
# 'DOF': 't',
# 'DIR': 'edgey',
# 'LOC': {'x': 999, 'y': 0, 'z': 0},
# 'VAL': [100.0],
# }
physicdata = {
"PHYSIC": {"DOMAIN": "thermal", # 'fluid' 'thermal'; "COUPLING": 'fsi'
"LOAD": [],
"BOUNDCOND": [bcNH1],
}
}
fea.Physic(physicdata)
previewset = {'RENDER': {'filename': 'heat_transfer', 'show': True, 'scale': 5, 'savepng': True, 'lines': True,
# 'plottags': {'line': True}
},
# 'LABELS': {'show': True, 'lines': True, 'scale': 1},
}
fea.PreviewAnalysis(previewset)
# #-------------------------------- SOLVER -------------------------------------#
solverset = {"STEPSET": {'type': 'table', # mode, freq, time ...
'start': 0,
'end': 1,
'step': 1},
'SYMM':True,
# 'MP':True,
}
solverdata = fea.Solve(solverset)
# print(solverdata['solution']['U'])
postprocset = {"SOLVERDATA": solverdata,
"COMPUTER": {'thermal': {'temp': True, 'heatflux': True}},
"PLOTSET": {'show': True, 'filename': 'SIM_TIRA_BIMETALICA_HEAT', 'savepng': True},
# "TRACKER": {'point': {'x': 0, 'y': 0, 'z': 0, 'dof':1}},
"OUTPUT": {'log': True, 'get': {'nelem': True, 'nnode': True}},
}
postprocdata = fea.PostProcess(postprocset)
### sova simulação, aplicando os valores do campo de calor da simulação anterior para gerar expansão termica no sólido
modeldata = {
"MESH": {
'TYPE': 'gmsh',
'filename': 'test_line_force_x',
# "meshimport": 'object_dir',
'pointlist': points,
'linelist': lines,
'planelist': plane,
# 'arc': arcs,
'meshconfig': {
'mesh': 'hexa8', #quad4
'sizeelement': 4,
'extrude': 20,
'meshmap': {'on': True,
'edge': 'all', #'all'
"numbernodes": [10],
}
}
},
"ELEMENT": {
'TYPE': 'structsolid',
'SHAPE': 'hexa8',
# 'INTGAUSS': 7,
},
"MATERIAL": {
"MAT": 'solidelastic',
"TYPE": 'isotropic',
"PROPMAT": [copper, steel],
},
"GEOMETRY": {
"GEO": 'thickness',
"PROPGEO": [geo, geo],
},
}
fea.Model(modeldata)
# fc = {
# 'TYPE': 'forcenode',
# 'DOF': 'fy',
# 'DIR': 'node',
# 'LOC': {'x': 200, 'y': LY, 'z': 0},
# 'VAL': [0.0],
# }
fc = {
'TYPE': 'forcenode',
'DOF': 'fy',
'DIR': 'edgex',
'LOC': {'x': 200, 'y': 999, 'z': 0},
'VAL': [0.0],
}
# bcb = {
# 'TYPE': 'fixed',
# 'DOF': 'uy',
# 'DIR': 'edgey',
# 'LOC': {'x': 999, 'y': 0, 'z': 0},
# }
bcl = {
'TYPE': 'fixed',
'DOF': 'full',
'DIR': 'surfyz',
'LOC': {'x': 0, 'y': 999, 'z': 999},
}
# bcr = {
# 'TYPE': 'fixed',
# 'DOF': 'full',
# 'DIR': 'edgex',
# 'LOC': {'x': LX, 'y': 999, 'z': 0},
# }
physicdata = {
"PHYSIC": {"DOMAIN": "structural",
"LOAD": [],
"BOUNDCOND": [bcl],
},
"COUPLING": {"TYPE": 'thermalstress', # fsi asi
"POST": [postprocdata]
},
}
fea.Physic(physicdata)
previewset = {'RENDER': {'filename': 'structural', 'show': True, 'scale': 10, 'savepng': True, 'lines': False,
# 'plottags': {'line': True}
},
# 'LABELS': {'show': True, 'lines': True, 'scale': 1},
}
fea.PreviewAnalysis(previewset)
# #-------------------------------- SOLVER -------------------------------------#
solverset = {"STEPSET": {'type': 'table', # mode, freq, time ...
'start': 0,
'end': 1,
'step': 1},
'SYMM':True,
# 'MP':True,
}
solverdata = fea.Solve(solverset)
postprocset = {"SOLVERDATA": solverdata,
"COMPUTER": {'structural': {'displ': True, 'stress': True}},
"PLOTSET": {'show': True, 'filename': 'SIM_TIRA_BIMETALICA_DISP', 'savepng': True},
# "TRACKER": {'point': {'x': 0, 'y': 0, 'z': 0, 'dof':1}},
"OUTPUT": {'log': True, 'get': {'nelem': True, 'nnode': True}},
}
postprocdata = fea.PostProcess(postprocset)
SET NEW MATERIAL: ORTHOTROPIC ELASTIC¶
SET NEW MATERIAL ORTHOTROPIC ELASTIC¶
import numpy as np
from myfempy import newAnalysis
from myfempy import SteadyStateLinear, SteadyStateLinearIterative
from myfempy import Material, PlaneStress
# ===============================================================================
# FEA
# ===============================================================================
fea = newAnalysis(SteadyStateLinear)
# #-------------------------------- PRE -------------------------------------#
mat = {
"NAME": "fibra_vidro", # Fonte: CALLISTER
"EXX": 45E3, # N/mm^2 --> MPa
"EYY": 10E3, # N/mm^2 --> MPa
"VXY": 0.3,
"VYZ": 0.4,
}
geo = {
"NAME": "espessura",
"THICKN": 5, # mm
}
# ===============================================================================
# SET NEW MATERIAL <ORTHOTROPIC ELASTIC>
# ===============================================================================
class UserNewMaterial(PlaneStress):
def __init__(self):
super().__init__()
def getMaterialSet():
matset = {
"mat": "planestress",
"type": "orthotropic",
}
return matset
def getElasticTensor(Model, element_number):
# # material elasticity
EXX = Model.tabmat[int(Model.inci[element_number, 2]) - 1]["EXX"]
EYY = Model.tabmat[int(Model.inci[element_number, 2]) - 1]["EYY"]
# material poisson ratio
VXY = Model.tabmat[int(Model.inci[element_number, 2]) - 1][ "VXY"]
VYZ = Model.tabmat[int(Model.inci[element_number, 2]) - 1][ "VYZ"]
# EXX = 45E3 # N/mm^2 --> MPa
# EYY = 12E3 # N/mm^2 --> MPa
# VXY = 0.23
# VYX = 0.66
S00 = EXX/(1-VXY*VYZ)
S01 = (VYZ*EXX)/(1-VXY*VYZ)
S10 = (VXY*EYY)/(1-VXY*VYZ)
S11 = EYY/(1-VXY*VYZ)
S22 = (EXX*EYY)/(EXX+EYY+2.0*EYY*VXY)
D = np.zeros((3, 3), dtype=np.float64)
D[0, 0] = S00
D[0, 1] = S01
D[1, 0] = S10
D[1, 1] = S11
D[2, 2] = S22
return D
# def getFailureCriteria(stress):
# stress
# return super().getFailureCriteria()
# ===============================================================================
# config model geometry --> mesh(gmsh)
points = [
[0, 0, 0],
[200, 0, 0],
[0, 10, 0],
[200, 10, 0],
[0, 20, 0],
[200, 20, 0]
]
lines = [[1, 2],
[2, 4],
[3, 4],
[3, 1],
[6, 4],
[5, 6],
[5, 3],
]
plane = [[1, 2, 3, 4], [3, 5, 6, 7]]
# #-------------------------------- THERMAL SIM -------------------------------------#
modeldata = dict()
modeldata["MESH"] = {
'TYPE': 'gmsh',
'filename': 'test_line_force_x',
'pointlist': points,
'linelist': lines,
'planelist': plane,
# 'arc': arcs,
'meshconfig': {
'mesh': 'quad4', #quad4
'sizeelement': 1,
'meshmap': {'on': True,
'edge': [[1,3,6], [2,4,5,7]], #'all'
"numbernodes": [20, 10],
}
}
}
modeldata["ELEMENT"] = {
'TYPE': 'structplane',
'SHAPE': 'quad4',
'INTGAUSS': 4,
}
modeldata["MATERIAL"] = {
"MAT": 'planestress',
"TYPE": 'usermaterial',
"PROPMAT": [mat, mat],
"CLASS": UserNewMaterial,
}
modeldata["GEOMETRY"] = {
"GEO": 'thickness',
"PROPGEO": [geo, geo],
}
fea.Model(modeldata)
f1 = {
'TYPE': 'forcenode',
'DOF': 'fx',
'DIR': 'edgex',
'LOC': {'x': 200, 'y': 999, 'z': 0},
'VAL': [10.0],
}
bc1 = {
'TYPE': 'fixed',
'DOF': 'full',
'DIR': 'edgex',
'LOC': {'x': 0, 'y': 999, 'z': 0},
}
physicdata = {
"PHYSIC": {"DOMAIN": "structural",
"LOAD": [f1],
"BOUNDCOND": [bc1],
},
}
fea.Physic(physicdata)
previewset = {'RENDER': {'filename': 'model_preview', 'show': True, 'scale': 2, 'savepng': True, 'lines': False,
'plottags': {'line': True}
},
# 'LABELS': {'show': True, 'lines': True, 'scale': 1},
}
fea.PreviewAnalysis(previewset)
# sys.exit()
# #-------------------------------- SOLVER -------------------------------------#
solverset = {"STEPSET": {'type': 'table',
'start': 0,
'end': 1,
'step': 1},
'SYMM':False,
}
solverdata = fea.Solve(solverset)
# print(solverdata['solution']['U'])
postprocset = {"SOLVERDATA": solverdata,
"COMPUTER": {'structural': {'displ': True, 'stress': True}},
"PLOTSET": {'show': True, 'filename': 'solution', 'savepng': True},
"OUTPUT": {'log': True, 'get': {'nelem': True, 'nnode': True}},
}
postprocdata = fea.PostProcess(postprocset)
Appendix¶
Useful Links¶
GMSH FreeCAD Python NumPy SciPy Computer-aided design Finite element method Young's modulus Poisson's ratio Shear modulus Density List of moments of inertia Boundary value problem Gaussian quadrature International System of Units
Axis Diretions¶
# |
# [Y]
# | P1 -- principal plane
# | P2 -- secondary plane
# |__edgey__
# /| |
# / | P1 |
# / | surfxy edgex
# / f| |
# / r |_________|____________[X]__
# | u / /
# |s z/ P2 /
# | y/ surfzx edgez
# | / /
# |/__________/
# /
# [Z]
#/
Cross Section Dimensions¶
Geometric Cross-Section Libraries
# :
# [Y]
# :
# ___ ___________:___________
# | |_______ : _______|
# | | : |
# | -->| : |<-------------(t)
# | | : |
# | | : |
# | | : |
# (h) | (CG).|..........[Z]..
# | | |
# | | |
# | | |
# | | |
# | _______| |_______ ___
# _|_ |_______________________| _|_(d)
# |-----------(b)---------|
Table 1 - Elements List¶
| element | id | description |
|---|---|---|
| block | 11 | Spring + Mass 1D-space 1-node_dofs |
| structbeam | 16 | Beam Structural Element 1D-space 6-node_dofs |
| structplane | 22 | Plane Structural Element 2D-space 2-node_dofs |
| structsolid | 33 | Solid Structural Element 3D-space 3-node_dofs |
| heatplane | 21 | Plane Heat Element 2D-space 1-node_dof |
| heatsolid | 31 | Solid Heat Element 3D-space 1-node_dofs |
Table 2 - Shape (Mesh) List¶
| mesh | id | description | Gauss Points Opt. |
|---|---|---|---|
| line2 | 21 | Line 2-Node Shape 2-nodes_conec 1-interpol_order | 1/ 2/ 4/ 8 |
| line3 | 32 | Line 3-Node Shape 3-nodes_conec 2-interpol_order | 1/ 2/ 4/ 8 |
| tria3 | 31 | Triangular 3-Node Shape 3-nodes_conec 1-interpol_order | 1/ 3/ 7 |
| tria6 | 62 | Triangular 6-Node Shape 6-nodes_conec 2-interpol_order | 1/ 3/ 7 |
| quad4 | 41 | Quadrilateral 4-Node Shape 4-nodes_conec 1-interpol_orde | 1/ 2/ 3/ 4/ 8/ 9 |
| quad8 | 82 | Quadrilateral 8-Node Shape 8-nodes_conec 2-interpol_order | 1/ 2/ 3/ 4/ 8/ 9 |
| hexa8 | 81 | Hexaedron 8-Node Shape 8-nodes_conec 1-interpol_order | 1/ 2/ 3/ 4/ 8/ 9 |
| tetr4 | 41 | Tetrahedron 4-Node Shape 4-nodes_conec 1-interpol_orde | 1/ 3/ 5 |
Table 3 - Solvers List¶
| solver | id | description | mandatory parameters | optional parameters |
|---|---|---|---|---|
| SteadyStateLinear | Steady State Linear Solver Class (Static Linear) | Model, Assembly: [stiffness, loads], ConstrainsDof: [freedof, constdof] | Assembly: [bcdirnh], Solverset: [nsteps] | |
| SteadyStateLinearIterative | Steady State Linear Iterative Solver Class (Static Linear) | Model, Assembly: [stiffness, loads], ConstrainsDof: [freedof, constdof] | Assembly: [bcdirnh], Solverset: [nsteps, tol, maxiter] | |
| DynamicEigenLinear | Dynamic Eigen (modal problem) Linear Solver Class | Model, Assembly: [stiffness, mass], ConstrainsDof: [freedof] | Solverset: [modeEnd:nsteps] | |
| DynamicHarmonicResponseLinear | Dynamic Harmonic Response Forced System Steady State Linear Solver Clas | Model, Assembly: [stiffness, mass, loads], ConstrainsDof: [freedof] | Solverset: [freqStart:start, freqEnd:end,freqStep:nsteps] | |
| [ dev ] TransientLinear | ||||
| [ adv ] StaticLinearCyclicSymmPlane | Static Linear Cyclic Symmetry Plane Solver Class | Model, Physic, Assembly: [stiffness, loads], ConstrainsDof: [freedof, fixedof(leftdof, rightdof, interdof)*] | Solverset: [nsteps] | |
| [ adv ] HomogenizationPlane | Homogenization Plane Solver Class | Model, Assembly: [stiffness, loads], ConstrainsDof: [freedof, fixedof(strain)*] | ||
| [ adv ] PhononicCrystalInPlane | Phononic Crystal In-Plane Solver Class | Model, Assembly: [stiffness, mass], ConstrainsDof: [freedof, constdof(bloch)*] | Solverset: [nsteps, IBZ*] | |
| [ dev ] TopologyOptimization | ||||
*Note: See Solver Set
Legends
- [ Model ]: Model object of analysis
- [ Physic ]: Physic object of analysis
- [ Assembly ]: Dictionary with the algebraic matrices
- [ stiffness ]: linear stiffness
- [ mass ]: mass matrix
- [ loads ]: vector of loads of the degrees of freedom
- [ bcdirnh ]: vector with Dirichlet boundary conditions (non-homogeneous)
- [ ConstrainsDof ]: dictionary with the problem's restrictions
- [ freedof ]: free d.o.f.
- [ fixeddof ]: fixed d.o.f.
- [ constdof ]: restricted d.o.f.
- [ Solverset ]: dictionary with commands for the solver
- [ nsteps ]: number of solver steps
- [ start ]: first step number
- [ end ]: last step number
- [ IBZ ]: Brillouin zone rad/s
Table 4 - Material List¶
| material | id | description |
|---|---|---|
| UniAxialStress | Uni-Axial Stress Isotropic Material Class | |
| PlaneStress | Plane Stress Isotropic Material Class | |
| PlaneStrain | Plane Strain Isotropic Material Class | |
| SolidElastic | Solid Stress Isotropic Material Class | |
| HeatPlane | Heat Plane Isotropic Material Class | |
| HeatSolid | Heat Solid Isotropic Material Class | |
| UserNewMaterial | User Defined Material API | |
Table 5 - Consistent Units¶
| Quantity | SI(m) | SI(mm) |
|---|---|---|
| length | m | mm |
| mass | kg | tonne |
| time | s | s |
| temperature | °C | °C |
| rotation | rad | rad |
| acceleration/ gravity | m/s^2 | mm/s^2 |
| density | kg/m^3 | tonne/mm^3 |
| force | N(kg⋅m⋅s^−2) | N |
| moment | N-m | N-mm |
| frequency | Hz(s^-1) = (rad/s)/(2.pi) | Hz |
| pressure/ stress | Pa(N/m^2) | MPa(N/mm^2) |
| energy | J(kg⋅m^2⋅s^−2) | mJ(tonne⋅mm^2⋅s^−2) |
| power | W(kg⋅m^2⋅s^-3) | mW(tonne⋅mm^2⋅s^-3) |
| thermal conductivity | W/(m.°C) | mW/(mm·°C) |
| thermal expansion | (m/m)/°C | (mm/mm)/°C |
| heat flux | W/m^2 | W/mm^2 |
| viscosity | Pa.s | MPa.s |
See International System of Units
Tag Legends¶
- [ adv ]: Inputs advanced options or external package
- [ dev ]: Inputs options in development (next update)
- [ old ]: Inputs of legacy/old version