# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
"""Samcef output file reader
"""
import subprocess
import os
import numpy as np
import shlex
import Muscat.Helpers.ParserHelper as PH
from Muscat.Helpers.Logger import Debug
from Muscat.IO.SamcefTools import PropertiesReader
import Muscat.MeshContainers.Mesh as UM
from Muscat.Types import MuscatFloat, MuscatIndex
import Muscat.MeshContainers.ElementsDescription as ED
samresExec = "samres"
[docs]
class SamcefOuputReader():
"""Samcef output Reader class
"""
def __init__(self):
super().__init__()
self.canHandleTemporal = True
self.subprocess = None
self.filename = ""
self.LCP = "me"
self.timeToRead = 0.
self.fieldsToRead = ["U"]
self.encoding = None
[docs]
def SetFileName(self, fileName):
"""Sets the name of the file to read
Parameters
----------
fileName : str
name of the file to read
"""
Debug(fileName)
if len(fileName) > 4 and fileName[-4:] == ".fac":
fileName = fileName[0:-4]
if len(fileName) > 3 and fileName[-3] == "_":
self.LCP = fileName[-2:]
fileName = fileName[0:-3]
self.filename = fileName
[docs]
def SetTimeToRead(self, time = None, timeIndex = None):
"""Sets the time at which the data is read
Parameters
----------
time : float, optional
time at which the data is read, by default None
timeIndex : int, optional
time index which the data is read, by default None
"""
if timeIndex is not None:
self.timeToRead = self.times[timeIndex]
elif time is not None:
self.timeToRead = PH.ReadFloat(time)
[docs]
def Open(self):
if self.subprocess is not None:
raise(Exception("File already open"))
casename = self.filename.split(".fac")[0]
args = [samresExec,"NOM="+casename, "LCP="+self.LCP]
Debug("Opening "+ str(casename))
self.subprocess = subprocess.Popen(args, bufsize=1, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE, universal_newlines=True)
#executable=None,stderr=None, preexec_fn=None, close_fds=True, shell=False, cwd=None, env=None, startupinfo=None, creationflags=0, restore_signals=True, start_new_session=False, pass_fds=(), *, encoding=None, errors=None)
[docs]
def OpenIfNeeded(self):
if self.subprocess == None:
self.Open()
return True
else:
return False
[docs]
def Close(self,needToClose=True):
if needToClose :
self.subprocess.terminate()
self.subprocess.communicate()
self.subprocess = None
[docs]
def GetAvailableTimes(self):
"""Returns available times to read
Returns
-------
list
time values available to read
"""
return self.times
[docs]
def ReadMeshStructure(self):
from Muscat.IO.SamcefReader import DatReader
dreader = DatReader()
mesh = dreader.Read(fileName=self.filename+".dat")
foid = dreader.filetointernalid
foid_element = dreader.filetointernalidElement
return mesh, foid, foid_element
[docs]
def Read(self):
"""Function that performs the reading of a samcef output file
Returns
-------
Mesh
output unstructured mesh object containing reading result
"""
Debug("start reading")
mesh, foid, foid_element = self.ReadMeshStructure()
Debug("reading mesh")
needClose = self.OpenIfNeeded()
name3 = ["X","Y","Z"]
name6 = ["XX","YY","ZZ","XY", "XZ", "YZ"]
names = {3:name3, 6:name6}
offsets= {}
cpt = 0
for elemtype, data in mesh.elements.items():
offsets[elemtype] = cpt
cpt += data.GetNumberOfElements()
for fname in self.fieldsToRead:
Debug("Reading "+ fname )
oid, ofield = self.ReadField(fname)
if InternalCodes[fname].on == "Nodes":
field = np.zeros( (mesh.GetNumberOfNodes(),ofield.shape[1]) )
for i in range(ofield.shape[0]):
for j in range(ofield.shape[1]):
field[foid[oid[i,j]],j ] = ofield[i,j]
storage = mesh.nodeFields
elif InternalCodes[fname].on == "Elements":
field = np.zeros( (mesh.GetNumberOfElements(),ofield.shape[1]) )
for i in range(ofield.shape[0]):
for j in range(ofield.shape[1]):
elem,cpt = foid_element[oid[i,j]]
enumber = offsets[elem.elementType]+cpt
field[enumber,j ] = ofield[i,j]
storage =mesh.elemFields
if field.shape[1] > 3:
for i in range(field.shape[1]):
storage[fname+names[field.shape[1]][i]] = field[:,i]
else:
storage[fname] = field
self.Close(needClose)
return mesh
[docs]
def ReadField(self, fieldname = None, time = None, timeIndex = None):
"""Function that performs the reading of a field in a samcef output file
Parameters
----------
fieldname : str, optional
name of the field to read, by default None
time : float, optional
time at which the data is read, by default None
timeIndex : int, optional
time index which the data is read, by default None
Returns
-------
np.ndarray, np.ndarray
indices and values of the field to read
"""
needToClose = self.OpenIfNeeded()
# documentation samres page 25
self.SetTimeToRead(time=time,timeIndex=timeIndex)
arg1,arg2,arg3,arg4,arg5,arg6 = InternalCodes[fieldname].GetArgs()
nbcomp= InternalCodes[fieldname].nbcomp
if self.timeToRead is not None:
arg5 = 'Time '+ str(self.timeToRead)
else:
arg5 = " "
cmd =r"$$GET_VALUE " + " ".join(['"'+str(x)+'"' for x in [arg1,arg2,arg3,arg4,arg5,arg6]])
self._writeCommand(cmd)
NRef = PH.ReadInt(self._readResponse())
refNames = []
for i in range(NRef):
refNames.append(self._readResponse())
SEL_SIZE = PH.ReadInt(self._readResponse())//nbcomp
entityname = self._readResponse()
oid = self._readNumpyArray((SEL_SIZE,nbcomp), dtype=int)
data = self._readNumpyArray((SEL_SIZE,nbcomp), dtype=float)
self.Close(needToClose)
if entityname not in ['"Node %"','"Element %"']:
raise(Exception("Dont know how to treat->"+ entityname+"<-"))
return oid, data
def _writeCommand(self,cmd):
Debug("Command:"+cmd)
self.subprocess.stdin.write(str(cmd)+os.linesep)
self.subprocess.stdin.flush()
self._skipResponse()
def _skipResponse(self,nblines=1):
Debug(f"Skiping {nblines} lines ")
for i in range(nblines):
self._readResponse()
def _readResponse(self):
data = self.subprocess.stdout.readline().rstrip()
if data.find("ERROR") >= 0:
raise(Exception(data))
return data
def _readNumpyArray(self, shape, dtype=float):
data = np.loadtxt(self.subprocess.stdout, dtype=dtype, max_rows=int(np.prod(shape)))
data.shape = shape
return data
from Muscat.IO.IOFactory import RegisterReaderClass
RegisterReaderClass(".fac",SamcefOuputReader)
[docs]
class SamcefData():
def __init__(self,arg1=" ",arg2=" ",arg3=" ",arg4=" ",name=" ",nbcomp=1):
self.arg1 = arg1
self.arg2 = arg2
self.arg3 = arg3
self.arg4 = arg4
self.nbcomp = nbcomp
[docs]
def GetArgs(self):
return (self.arg1,self.arg2,self.arg3,self.arg4," "," ")
[docs]
class SamcefDataN(SamcefData):
def __init__(self,*k,**nk):
super().__init__(*k,**nk)
self.on = "Nodes"
[docs]
class SamcefDataE(SamcefData):
def __init__(self,*k,**nk):
super().__init__(*k,**nk)
self.on = "Elements"
[docs]
class SamcefDataG(SamcefData):
def __init__(self,*k,**nk):
super().__init__(*k,**nk)
self.on = "Global"
InternalCodes = {}
InternalCodes["Time"] = SamcefData(arg1='Code 26',name="Time")
InternalCodes["U"] = SamcefDataN(arg1='Code 163',name="U",nbcomp=3)
InternalCodes["U0"] = SamcefDataN(arg1='Code 163',arg2="Component 1", name="U0")
InternalCodes["U1"] = SamcefDataN(arg1='Code 163',arg2="Component 2", name="U1")
InternalCodes["U2"] = SamcefDataN(arg1='Code 163',arg2="Component 3", name="U2")
InternalCodes["Reaction"] = SamcefDataN(arg1='Code 221',name="Reaction",nbcomp=3)
InternalCodes["Nodal_Stress"] = SamcefDataN(arg1='Code 1411',name="Nodal_Stress",nbcomp=6)
InternalCodes["Element_Stress"] = SamcefDataE(arg1='Code 3411',name="Element_Stress",nbcomp=6)
InternalCodes["Element_Stress_VM"] = SamcefDataE(arg1='Code 3411',arg2="Von Mises",name="Element_Stress_VM",nbcomp=1)
"""
"Sdb 4" "[Sdb] Element Mass" " " ".SUPPORT_0102"
"Sdb 3" "[Sdb] Element Volume" " " ".SUPPORT_0102"
"Sdb 2" "[Sdb] Element Surface" " " ".SUPPORT_0102"
"Sdb 1" "[Sdb] Element Length" " " ".SUPPORT_0102"
"Code 26" "[Code 26] Time" " " ".SUPPORT_0050"
"Code 163" "[Code 163] Nodal displacements (DX,DY,DZ)" ".TOSCALAR_0003" ".SUPPORT_0003"
"Code 191" "[Code 191] Nodal rotation (RX,RY,RZ)" ".TOSCALAR_0009" ".SUPPORT_0009"
"Code 221" "[Code 221] Nodal reaction" ".TOSCALAR_0010" ".SUPPORT_0010"
"Code 901" "[Code 901] Time" ".TOSCALAR_0901" ".SUPPORT_0901"
"Code 981" "[Code 981] Kinetic energy" ".TOSCALAR_0981" ".SUPPORT_0981"
"Code 982" "[Code 982] Potential energy" ".TOSCALAR_0982" ".SUPPORT_0982"
"Code 984" "[Code 984] External forces work" ".TOSCALAR_0984" ".SUPPORT_0984"
"Code 985" "[Code 985] Time steps" ".TOSCALAR_0985" ".SUPPORT_0985"
"Code 986" "[Code 986] Time integration errors" ".TOSCALAR_0986" ".SUPPORT_0986"
"Code 1305" "[Code 1305] Contact pressure" " " ".SUPPORT_0026"
"Code 1306" "[Code 1306] Friction stress" " " ".SUPPORT_0026"
"Code 1307" "[Code 1307] Normal distance" " " ".SUPPORT_0026"
"Code 1342" "[Code 1342] Sliding velocity direction 1" " " ".SUPPORT_0026"
"Code 1343" "[Code 1343] Sliding velocity direction 2" " " ".SUPPORT_0026"
"Code 1344" "[Code 1344] Sliding" " " ".SUPPORT_0026"
"Code 1411" "[Code 1411] Stress tensor" ".TOSCALAR_0013" ".SUPPORT_0013"
"Code 1421" "[Code 1421] Strain tensor" ".TOSCALAR_0013" ".SUPPORT_0013"
"Code 1431" "[Code 1431] 2D linear stress tensor" ".TOSCALAR_0015" ".SUPPORT_0015"
"Code 1437" "[Code 1437] Normal force" ".TOSCALAR_0015" ".SUPPORT_0015"
"Code 1445" "[Code 1445] 2D strain tensor (upper skin)" ".TOSCALAR_0015" ".SUPPORT_0015"
"Code 1446" "[Code 1446] 2D strain tensor (lower skin)" ".TOSCALAR_0015" ".SUPPORT_0015"
"Code 1524" "[Code 1524] Force vector" ".TOSCALAR_0025" ".SUPPORT_0025"
"Code 1525" "[Code 1525] Moment vector" ".TOSCALAR_0025" ".SUPPORT_0025"
"Code 3408" "[Code 3408] 2D Cauchy stress tensor (low. GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3409" "[Code 3409] 2D Biot stress tensor (low. GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3411" "[Code 3411] Stress tensor" ".TOSCALAR_0014" ".SUPPORT_0014"
"Code 3418" "[Code 3418] 2D Cauchy strain tensor (low. GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3419" "[Code 3419] 2D Biot strain tensor (low. GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3431" "[Code 3431] 2D linear stress tensor" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3448" "[Code 3448] 2D Cauchy stress tensor (middle GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3449" "[Code 3449] 2D Biot stress tensor (middle GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3452" "[Code 3452] Cauchy strain tensor" ".TOSCALAR_0014" ".SUPPORT_0014"
"Code 3478" "[Code 3478] 2D Cauchy stress tensor (upper GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3479" "[Code 3479] 2D Biot stress tensor (upper GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3488" "[Code 3488] 2D Cauchy strain tensor (upper GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 3489" "[Code 3489] 2D Biot strain tensor (upper GP)" ".TOSCALAR_0006" ".SUPPORT_0006"
"Code 9821" "[Code 9821] Resultants - 3 components" ".TOSCALAR_9821" ".SUPPORT_9821"
"Code 9822" "[Code 9822] Number of plastic elements" ".TOSCALAR_9822" ".SUPPORT_9822"
"Code 9823" "[Code 9823] Number of damaged elements" ".TOSCALAR_9823" ".SUPPORT_9823"
"Code 9824" "[Code 9824] Number of nodes in contact" ".TOSCALAR_9824" ".SUPPORT_9824"
"Code 9850" "[Code 9850] Number of Cycles" ".TOSCALAR_9850" ".SUPPORT_9850"
"""
[docs]
def ReadFieldFromDataBase(fieldname, field, time=None):
"""Function API for reading a field in a samcef output file
Parameters
----------
fieldname : str
name of the file to be read
field : str
name of the field to read in the file
time : float, optional
time at which the data is read, by default None
Returns
-------
np.ndarray, np.ndarray
indices and values of the field to read
"""
reader = SamcefOuputReader()
reader.SetFileName(fieldname)
reader.Open()
return reader.ReadField(field,time=time)
[docs]
class SamcefOutputBaconReader(PropertiesReader, SamcefOuputReader):
"""
Reads and parses mesh and element data from SAMCEF output files using baconpost.
This class provides an interface to process SAMCEF output data, extract the mesh
structure. This class inherits from PropertiesReader and can be used to extract
other relevant information (element attributes, mechanical boundary conditions).
Examples
--------
>>> reader = SamcefOutputBaconReader(folderpath=tempdir, projectname=projectname)
>>> reader.SetFileName(filename)
>>> reader.fieldsToRead = ["U"]
>>> reader.LCP = "as"
>>> reader.timeToRead = None
>>> mesh = reader.Read()
>>> element_id_to_element_info = reader.ReadElementAttributes() # uses class PropertiesReader
>>> node_id_to_clm, element_id_to_clm, structure_id_to_clm = reader.ReadMechanicalBoundaryConditions() # uses class PropertiesReader
>>> group_id_to_infos = reader.SelRead() # uses class PropertiesReader
"""
def __init__(self, folderpath=None, projectname=None):
super().__init__(folderpath=folderpath, projectname=projectname)
[docs]
def SetFileName(self, fileName):
SamcefOuputReader.SetFileName(self, fileName)
[docs]
def GetInternalNodeNumberFromOriginalID(self, theList):
try:
return [self.filetointernalid[x] for x in theList]
except KeyError as e:
raise IndexError("ERROR!!! Node with number " + str(e) + " not in file")
[docs]
def decoupe_par_zero(self, lst):
result = []
current = []
for x in lst:
if x == 0:
if current: # si on a quelque chose avant le 0
result.append(current)
current = []
else:
current.append(x)
if current: # ne pas oublier la dernière séquence
result.append(current)
return result
[docs]
def ReadMeshStructure(self):
"""
Read a SAMCEF output file and construct a mesh object.
This function reads only the nodes and elements from a SAMCEF output file
using baconpost and constructs a mesh object. The function makes certain
assumptions about element types.
Returns
-------
Mesh
Output mesh object.
"""
#read nodes and elements
original_ids, xs, ys, zs = self.NodesRead()
element_id_to_nodes, original_ids_elements = self.ElementRead()
#constructing the muscat mesh
res = UM.Mesh()
#node id in samcef to node id in muscat dictionnary
self.filetointernalid = {original_ids[i]:i for i in range(len(original_ids))}
self.filetointernalidElement = {original_ids_elements[i]:i for i in range(len(original_ids_elements))}
for element_id in element_id_to_nodes:
connectivity_original_ids = element_id_to_nodes[element_id]
#negative sign for interface nodes
connectivity_original_ids = [abs(x) for x in connectivity_original_ids]
#connectivity_original_ids are the nodes composing the element element_id
p2 = self.decoupe_par_zero(connectivity_original_ids)
if len(p2) == 2:
connectivity_original_ids = p2[0] + p2[1]
if len(p2[0]) == 3 and len(p2[1]) == 1:
elements = res.GetElementsOfType(ED.Tetrahedron_4)
elif len(p2[0]) == 4 and len(p2[1]) == 4:
elements = res.GetElementsOfType(ED.Hexahedron_8)
elif len(p2[0]) == 3 and len(p2[1]) == 3:
elements = res.GetElementsOfType(ED.Wedge_6)
else:
raise (Exception("type of element no coded" + str(p2)))
elif len(p2) == 1:
connectivity_original_ids = p2[0]
if len(p2[0]) == 3:
elements = res.GetElementsOfType(ED.Triangle_3)
elif len(p2[0]) == 4:
elements = res.GetElementsOfType(ED.Quadrangle_4)
else:
raise (Exception("type of element no coded" + str(p2)))
else:
raise (Exception("type of element no coded" + str(p2)))
connectivity = self.GetInternalNodeNumberFromOriginalID(connectivity_original_ids)
elements.AddNewElement(connectivity, element_id)
res.SetNodes(np.array([xs, ys, zs], dtype=MuscatFloat).T)
res.originalIDNodes = np.array(original_ids, dtype=MuscatIndex)
res.PrepareForOutput()
foid = self.filetointernalid
foid_element = self.filetointernalidElement
return res, foid, foid_element
RegisterReaderClass(".sdb",SamcefOutputBaconReader)
[docs]
def CheckIntegrity():
try:
subprocess.run("samcef", shell=True, check=True)
except:
return "skip"
data = u""".ssls06-m020 &
!
! THIN CYLINDER UNDER NORMAL PRESSURE
! Linear static analysis
!
.DEL.*
ABRE '/L' '4' ! Length, m
ABRE '/R' '1' ! Radius, m
ABRE '/T' '0.02' ! Thickness, m
!
ABRE '/Young' '2.1E11' ! Young modulus
ABRE '/nu' '0.3' ! Poisson ratio
ABRE '/P' '10000' ! Pressure, Pa
ABRE '/NL' '4' ! Number of elements
!
! A. GEOMETRY
! ===========
! A.1. 3D Geometry
! ----------------
.3POINT I 1 X 0. Y 0. Z 0.
I 2 X 0. Y (/R) Z 0.
I 3 X 0. Y -(/R) Z 0.
.ARC I 1 CENTRE 1 RAYON (/R)
!
! A.2. Mesher
! -----------
.CONTOUR AUTO OUVERT
.DOMAINE AUTO
!
.GEN DEGRE 1
MAILLE 1
GENERE
.EXT ATT 1
TZ (/L:2.) ELEMENT (/NL)
EXECUTE 1
!
!
! B. Hypothesis and material
! ==========================
.HYP MINDLIN
.MAT I 1
Yt (/Young)
Nt (/nu)
.PHP THICK ATT 1 V (/T)
.AEL MAT 1
!
! C. Boundary conditions
! ======================
! C.1. Fixations
! --------------
!.AXL LIGNE 1
! AXE 1 NORMAL
! AXE 2 TANGENT
.CLM FIX LIGNE 1 COMPO 3 4 5
FIX NOEUD 84 COMPO 1 2 6
!
! C.2. internal pressure
! ----------------------
.CLM ATT 1 PRESS V (/P)
!
! D. Run management
! =================
! (standard)
!
! E. Output control
! =================
! (standard + example of code removal)
.sai archive struct stype -1251 -1438
EXIT
!
! POSTPROCESSING
! ==============
!
.POST &
.DEL.*
.DOC DB
.DES
GRAP VISEE 1 1 1
CODE 163 COMP 1 ! Horizontal displacement
TRI 2 LIST 3
AFF VI
WAIT
!
CODE 163 COMP 3 ! Vertical displacement
TRI 2 LIST 3
AFF VI
WAIT
!
CODE 1431 COMP 1 ! Stress
AFF VI
WAIT
COMP 2
AFF VI
RETURN
"""
from Muscat.Helpers.CheckTools import SkipTest
from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory
from Muscat.Helpers.IO.Which import Which
if SkipTest("SAMCEF_NO_FAIL"): return "ok"
samresExecName = Which(samresExec)
if samresExecName == '' or samresExecName == None:
return "skip"
tempdir = TemporaryDirectory.GetTempPath()+"test/"
os.makedirs(tempdir, exist_ok=True)
projectname = "test"
filename = tempdir+projectname
f = open(filename+".dat", "w")
f.write(data)
f.close()
#Launch Samcef simulation (creates a .sdb file)
cmd = (
f'cd {tempdir} && samcef ba,as test n 1 zone=1000000000 >> test.log <<INPUTCMD\n'
f';INPUT .ASEF\n'
f';.FIN 1\n'
f'INPUTCMD\n'
)
subprocess.run(cmd, shell=True, check=True)
reader = SamcefOutputBaconReader(folderpath=tempdir, projectname=projectname)
reader.SetFileName(filename)
reader.LCP = "as"
reader.timeToRead = None
mesh = reader.Read()
element_id_to_element_info = reader.ReadElementAttributes()
element_id_to_nodes, _ = reader.ElementRead()
node_id_to_clm, element_id_to_clm, structure_id_to_clm = reader.ReadMechanicalBoundaryConditions()
group_id_to_infos = reader.SelRead()
return "ok"
if __name__ == '__main__':
print(CheckIntegrity())# pragma: no cover