Read Calculix Result files (*.frd) mesh

This commit is contained in:
jriegel 2013-12-30 22:33:05 +01:00
parent 360b190f03
commit c1a3b1434a
3 changed files with 192 additions and 62 deletions

View File

@ -64,7 +64,8 @@ SOURCE_GROUP("Module" FILES ${Mod_SRCS})
SET(FemScripts_SRCS
convert2TetGen.py
FemLib.py
FemLib.py
CalculixLib.py
MechanicalAnalysis.py
MechanicalMaterial.py
)

View File

@ -181,21 +181,42 @@ PyObject* FemMeshPy::compute(PyObject *args)
PyObject* FemMeshPy::addNode(PyObject *args)
{
double x,y,z;
if (!PyArg_ParseTuple(args, "ddd",&x,&y,&z))
return 0;
int i = -1;
if (PyArg_ParseTuple(args, "ddd",&x,&y,&z)){
try {
SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
SMDS_MeshNode* node = meshDS->AddNode(x,y,z);
if (!node)
throw std::runtime_error("Failed to add node");
return Py::new_reference_to(Py::Int(node->GetID()));
}
catch (const std::exception& e) {
PyErr_SetString(PyExc_Exception, e.what());
return 0;
}
}
PyErr_Clear();
try {
SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
SMDS_MeshNode* node = meshDS->AddNode(x,y,z);
if (!node)
throw std::runtime_error("Failed to add node");
return Py::new_reference_to(Py::Int(node->GetID()));
}
catch (const std::exception& e) {
PyErr_SetString(PyExc_Exception, e.what());
return 0;
if (PyArg_ParseTuple(args, "dddi",&x,&y,&z,&i)){
try {
SMESH_Mesh* mesh = getFemMeshPtr()->getSMesh();
SMESHDS_Mesh* meshDS = mesh->GetMeshDS();
SMDS_MeshNode* node = meshDS->AddNodeWithID(x,y,z,i);
if (!node)
throw std::runtime_error("Failed to add node");
return Py::new_reference_to(Py::Int(node->GetID()));
}
catch (const std::exception& e) {
PyErr_SetString(PyExc_Exception, e.what());
return 0;
}
}
PyErr_SetString(PyExc_TypeError, "addNode() accepts:\n"
"-- addNode(x,y,z)\n"
"-- addNode(x,y,z,ElemId)\n");
return 0;
}
PyObject* FemMeshPy::addEdge(PyObject *args)
@ -356,24 +377,47 @@ PyObject* FemMeshPy::addVolume(PyObject *args)
}
SMDS_MeshVolume* vol=0;
switch(Nodes.size()){
case 4:
vol = meshDS->AddVolume(Nodes[0],Nodes[1],Nodes[2],Nodes[3]);
if (!vol)
throw std::runtime_error("Failed to add Tet4 volume");
break;
case 8:
vol = meshDS->AddVolume(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7]);
if (!vol)
throw std::runtime_error("Failed to add Tet10 volume");
break;
case 10:
vol = meshDS->AddVolume(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7],Nodes[8],Nodes[9]);
if (!vol)
throw std::runtime_error("Failed to add Tet10 volume");
break;
if(ElementId != -1) {
switch(Nodes.size()){
case 4:
vol = meshDS->AddVolumeWithID(Nodes[0],Nodes[1],Nodes[2],Nodes[3],ElementId);
if (!vol)
throw std::runtime_error("Failed to add Tet4 volume");
break;
case 8:
vol = meshDS->AddVolumeWithID(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7],ElementId);
if (!vol)
throw std::runtime_error("Failed to add Tet10 volume");
break;
case 10:
vol = meshDS->AddVolumeWithID(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7],Nodes[8],Nodes[9],ElementId);
if (!vol)
throw std::runtime_error("Failed to add Tet10 volume");
break;
default: throw std::runtime_error("Unknown node count, [4|5|6|8|10|13|18] are allowed"); //unknown face type
}
}else{
switch(Nodes.size()){
case 4:
vol = meshDS->AddVolume(Nodes[0],Nodes[1],Nodes[2],Nodes[3]);
if (!vol)
throw std::runtime_error("Failed to add Tet4 volume");
break;
case 8:
vol = meshDS->AddVolume(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7]);
if (!vol)
throw std::runtime_error("Failed to add Tet10 volume");
break;
case 10:
vol = meshDS->AddVolume(Nodes[0],Nodes[1],Nodes[2],Nodes[3],Nodes[4],Nodes[5],Nodes[6],Nodes[7],Nodes[8],Nodes[9]);
if (!vol)
throw std::runtime_error("Failed to add Tet10 volume");
break;
default: throw std::runtime_error("Unknown node count, [4|5|6|8|10|13|18] are allowed"); //unknown face type
}
default: throw std::runtime_error("Unknown node count, [4|5|6|8|10|13|18] are allowed"); //unknown face type
}
return Py::new_reference_to(Py::Int(vol->GetID()));

View File

@ -1,5 +1,6 @@
#***************************************************************************
#* *
#* Copyright (c) 2013 - Joachim Zettler *
#* Copyright (c) 2013 - Juergen Riegel <FreeCAD@juergen-riegel.net> *
#* *
#* This program is free software; you can redistribute it and/or modify *
@ -21,42 +22,126 @@
#***************************************************************************
import FreeCAD,os
__title__="FreeCAD Calculix library"
__author__ = "Juergen Riegel "
__url__ = "http://www.freecadweb.org"
if open.__module__ == '__builtin__':
pyopen = open # because we'll redefine open below
# read a calculix result file and extract the nodes, displacement vectores and stress values.
def readResult(frd_input) :
input = open(frd_input,"r")
nodes_x = []
nodes_y = []
nodes_z = []
disp_x = []
disp_y = []
disp_z = []
displaced_nodes_x = []
displaced_nodes_y = []
displaced_nodes_z = []
input = pyopen(frd_input,"r")
nodes = {}
disp = {}
stress = {}
elements = {}
disp_found = False
nodes_found = True
nodes_found = False
stress_found = False
elements_found = False
elem = -1
while True:
line=input.readline()
if not line: break
#first lets extract the node and coordinate information from the results file
if nodes_found and (line[1:3] == "-1"):
nodes_x.append(float(line[13:25]))
nodes_y.append(float(line[25:37]))
nodes_z.append(float(line[37:49]))
#Check if we found displacement section
if line[5:9] == "DISP":
disp_found = True
#we found a displacement line in the frd file
if disp_found and (line[1:3] == "-1"):
disp_x.append(float(line[13:25]))
disp_y.append(float(line[25:37]))
disp_z.append(float(line[37:49]))
#Check for the end of a section
if line[1:3] == "-3":
#the section with the displacements and the nodes ended
disp_found = False
nodes_found = False
line=input.readline()
if not line: break
#Check if we found nodes section
if line[4:6] == "2C":
nodes_found = True
#first lets extract the node and coordinate information from the results file
if nodes_found and (line[1:3] == "-1"):
elem = int(line[4:13])
nodes_x = float(line[13:25])
nodes_y = float(line[25:37])
nodes_z = float(line[37:49])
nodes[elem] = FreeCAD.Vector(nodes_x,nodes_y,nodes_z)
#Check if we found nodes section
if line[4:6] == "3C":
elements_found = True
#first lets extract element number
if elements_found and (line[1:3] == "-1"):
elem = int(line[4:13])
#then the 10 id's for the Tet10 element
if elements_found and (line[1:3] == "-2"):
node_id_2 = int(line[3:13])
node_id_1 = int(line[13:23])
node_id_3 = int(line[23:33])
node_id_4 = int(line[33:43])
node_id_5 = int(line[43:53])
node_id_7 = int(line[53:63])
node_id_6 = int(line[63:73])
node_id_9 = int(line[73:83])
node_id_8 = int(line[83:93])
node_id_10 = int(line[93:103])
elements[elem] = (node_id_1,node_id_2,node_id_3,node_id_4,node_id_5,node_id_6,node_id_7,node_id_8,node_id_9,node_id_10)
#Check if we found displacement section
if line[5:9] == "DISP":
disp_found = True
#we found a displacement line in the frd file
if disp_found and (line[1:3] == "-1"):
elem = int(line[4:13])
disp_x = float(line[13:25])
disp_y = float(line[25:37])
disp_z = float(line[37:49])
disp[elem] = FreeCAD.Vector(disp_x,disp_y,disp_z)
if line[5:11] == "STRESS":
stress_found = True
#we found a displacement line in the frd file
if stress_found and (line[1:3] == "-1"):
elem = int(line[4:13])
stress_1 = float(line[13:25])
stress_2 = float(line[25:37])
stress_3 = float(line[37:49])
stress_4 = float(line[49:61])
stress_5 = float(line[61:73])
stress_6 = float(line[73:85])
stress[elem] = (stress_1,stress_2,stress_3,stress_4,stress_5,stress_6)
#Check for the end of a section
if line[1:3] == "-3":
#the section with the displacements and the nodes ended
disp_found = False
nodes_found = False
stress_found = False
elements_found = False
input.close()
FreeCAD.Console.PrintLog('Read Calculix result: ' + `len(nodes)` + ' Nodes, ' + `len(disp)` + ' Displacements and ' + `len(stress)` + ' Stress values\n')
return {'Nodes':nodes,'Tet10Elem':elements,'Displacement':disp,'Stress':stress}
def importFrd(filename):
m = readResult(filename);
if(len(m) > 0):
import Fem
if(m.has_key('Tet10Elem') and m.has_key('Nodes') ):
mesh = Fem.FemMesh()
nds = m['Nodes']
for i in nds:
n = nds[i]
mesh.addNode(n[0],n[1],n[2],i)
elms = m['Tet10Elem']
for i in elms:
e = elms[i]
mesh.addVolume([e[0],e[1],e[2],e[3],e[4],e[5],e[6],e[7],e[8],e[9]],i)
Fem.show(mesh)
def insert(filename,docname):
"called when freecad wants to import a file"
try:
doc = FreeCAD.getDocument(docname)
except:
doc = FreeCAD.newDocument(docname)
FreeCAD.ActiveDocument = doc
importFrd(filename)
def open(filename):
"called when freecad opens a file"
docname = os.path.splitext(os.path.basename(filename))[0]
insert(filename,docname)