Source code for plugins.ccxdat

#
##
##  SPDX-FileCopyrightText: © 2007-2023 Benedict Verhegghe <bverheg@gmail.com>
##  SPDX-License-Identifier: GPL-3.0-or-later
##
##  This file is part of pyFormex 3.4  (Thu Nov 16 18:07:39 CET 2023)
##  pyFormex is a tool for generating, manipulating and transforming 3D
##  geometrical models by sequences of mathematical operations.
##  Home page: https://pyformex.org
##  Project page: https://savannah.nongnu.org/projects/pyformex/
##  Development: https://gitlab.com/bverheg/pyformex
##  Distributed under the GNU General Public License version 3 or later.
##
##  This program is free software: you can redistribute it and/or modify
##  it under the terms of the GNU General Public License as published by
##  the Free Software Foundation, either version 3 of the License, or
##  (at your option) any later version.
##
##  This program is distributed in the hope that it will be useful,
##  but WITHOUT ANY WARRANTY; without even the implied warranty of
##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##  GNU General Public License for more details.
##
##  You should have received a copy of the GNU General Public License
##  along with this program.  If not, see http://www.gnu.org/licenses/.
##

"""
Interface with Calculix FE result files (.dat).

"""
import re

import numpy as np

import pyformex.arraytools as at

from pyformex.plugins.fe_post import FeResult
from pyformex.mesh import Mesh


_re_float = re.compile(r"[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)?")
_re_header = re.compile(
    r" *(?P<name>[^ ]+) *\(.*\) +for set (?P<set>[^ ]+) *and time *(?P<time>[^ ]+)")


######################### functions #############################

def skipBlankLine(fil):
    line = next(fil).strip()
    if len(line) > 0:
        raise ValueError("Expected a blank line")


# TODO: this could be speeded up as done in tetgen.py
[docs]def readDispl(fil, nnodes, nres): """Read displacements from a Calculix .dat file""" values = np.zeros((nnodes, nres), dtype=at.Float) for line in fil: if len(line.strip()) == 0: break s = line.split() i = int(s[0]) x = [float(f) for f in s[1:]] values[i-1] = x return values
[docs]def readStress(fil, nelems, ngp, nres): """Read stresses from a Calculix .dat file""" values = np.zeros((nelems, ngp, nres), dtype=at.Float) for line in fil: if len(line.strip()) == 0: break s = line.split() i, j = [int(v) for v in s[:2]] x = [float(v) for v in s[2:]] values[i-1, j-1] = x return values
[docs]def readResults(fn, DB, nnodes, nelems, ngp): """Read Calculix results file for nnodes, nelems, ngp Add results to the specified DB """ fil = open(fn, 'r') result = {} step = 1 for line in fil: m = _re_header.match(line) if m: m = m.groupdict() name = m['name'][:5] setname = m['set'] time = m['time'] print("Match %s %s %s" % (name, setname, time)) skipBlankLine(fil) if name == 'displ': result[name] = readDispl(fil, nnodes, 3) elif name == 'stres': result[name] = readStress(fil, nelems, ngp, 6) if 'displ' in result and 'stres' in result: addFeResult(DB, step, time, result) result = {} step += 1 return result
[docs]def createResultDB(model): """Create a results database for the given FE model""" DB = FeResult() DB.nodes = model.coords DB.nnodes = model.coords.shape[0] DB.nodid = np.arange(DB.nnodes) DB.elems = dict(enumerate(model.elems)) DB.nelems = model.celems[-1] DB.Finalize() return DB
[docs]def addFeResult(DB, step, time, result): """Add an FeResult for a time step to the result DB This is currently 2D only """ print("Storing result for step %s, time %s" % (step, time)) DB.Increment(step, 0) DB.R['TIME'] = time if 'displ' in result: DB.datasize['U'] = result['displ'].shape[1] DB.R['U'] = result['displ'] if 'stres' in result: try: stress = result['stres'] # CALCULIX HAS NO 2D: keep only half of GP's ngp = stress.shape[1]//2 stress = stress[:, :ngp, :] print("Reduced stresses: %s" % str(stress.shape)) mesh = Mesh(DB.nodes, DB.elems[0], eltype='quad4') gprule = [2, 2] stress = computeAveragedNodalStresses(mesh, stress, gprule) DB.datasize['S'] = result['stres'].shape[1] DB.R['S'] = stress except Exception: print("Error importing stresses") return DB
[docs]def computeAveragedNodalStresses(M, data, gprule): """Compute averaged nodal stresses from GP stresses in 2D quad8 mesh""" ############################ # Load the needed calpy modules from pyformex.plugins import calpy_itf calpy_itf.check() gprule = (2, 2) # integration rule: minimum (1,1), maximum (5,5) Q = calpy_itf.QuadInterpolator(M.nelems(), M.nplex(), gprule) ngp = np.prod(gprule) # number of datapoints per element print("Number of data points per element: %s" % ngp) print("Original element data: %s" % str(data.shape)) # compute the data at the nodes, per element endata = Q.GP2Nodes(data) print("Element nodal data: %s" % str(endata.shape)) # compute nodal averages nodata = Q.NodalAvg(M.elems+1, endata, M.nnodes()) print("Average nodal data: %s" % str(nodata.shape)) # extract the colors per element return nodata
# End