diff --git a/src/ifc2ca/Dockerfile b/src/ifc2ca/Dockerfile new file mode 100644 index 0000000000..8b0256b602 --- /dev/null +++ b/src/ifc2ca/Dockerfile @@ -0,0 +1,10 @@ +# Use Miniconda base image +FROM continuumio/miniconda3:4.10.3 + +# Update Conda, install necessary libraries, and then install Mamba +RUN conda update -n base -c defaults conda && \ + conda install libarchive -c conda-forge -y && \ + conda install mamba -c conda-forge -y + +# Install Code_Aster and Python dependencies with Conda +RUN mamba install -c conda-forge code-aster python=3.10 -y diff --git a/src/ifc2ca/__init__.py b/src/ifc2ca/__init__.py new file mode 100644 index 0000000000..8dd03376e1 --- /dev/null +++ b/src/ifc2ca/__init__.py @@ -0,0 +1,19 @@ +# Ifc2CA - IFC Code_Aster utility +# Copyright (C) 2020, 2021, 2023, 2024 Ioannis P. Christovasilis +# +# This file is part of Ifc2CA. +# +# Ifc2CA is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Ifc2CA 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Ifc2CA. If not, see . + +from .ifc2ca import Ifc2CA diff --git a/src/ifc2ca/_deprecated/ca2ifc.py b/src/ifc2ca/_deprecated/ca2ifc.py new file mode 100644 index 0000000000..7847a88e80 --- /dev/null +++ b/src/ifc2ca/_deprecated/ca2ifc.py @@ -0,0 +1,525 @@ + +# Ifc2CA - IFC Code_Aster utility +# Copyright (C) 2020, 2021 Ioannis P. Christovasilis +# +# This file is part of Ifc2CA. +# +# Ifc2CA is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Ifc2CA 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Ifc2CA. If not, see . + +import json +import ifcopenshell +import os +from datetime import datetime + + +class CA2IFC: + def __init__(self, inputFilename, outputFilename): + self.inputFilename = inputFilename + self.outputFilename = outputFilename + self.data = None + self.f = None + self.reps = {} + self.origin = None + self.xAxis = None + self.yAxis = None + self.zAxis = None + + def convert(self): + # load json file + with open(self.inputFilename) as dataFile: + self.data = json.load(dataFile) + + # initiate ifc file + self.f = ifcopenshell.file() + + # create header + self.create_header() + + # create global axes + globalAxes = self.create_global_axes() + localPlacement = self.f.createIfcLocalPlacement(None, globalAxes) + + # TODO: create units + lengthUnit = self.f.createIfcSIUnit(None, "LENGTHUNIT", None, "METRE") + unitAssignment = self.f.createIfcUnitAssignment((lengthUnit,)) + + # create owner history + ownerHistory = self.create_owner_history() + + # create representations and subrepresentations + self.reps = self.create_reference_subrep(globalAxes) + + # create project and model + project = self.f.createIfcProject( + self.guid(), ownerHistory, "A Project", None, None, None, None, (self.reps["model"],), unitAssignment + ) + model = self.f.createIfcStructuralAnalysisModel( + self.guid(), + ownerHistory, + self.data["name"], + None, + None, + "NOTDEFINED", + globalAxes, + None, + None, + localPlacement, + ) + self.f.createIfcRelDeclares(self.guid(), ownerHistory, None, None, project, (model,)) + + # create materials + ifcMaterials = [None for _ in range(len(self.data["db"]["materials"]))] + for i, material in enumerate(self.data["db"]["materials"]): + ifcMaterials[i] = self.create_material(material) + + # create profiles + ifcProfiles = [None for _ in range(len(self.data["db"]["profiles"]))] + for i, profile in enumerate(self.data["db"]["profiles"]): + ifcProfiles[i] = self.create_profile(profile) + + # create material-profile sets + mpSets = list( + set([el["material"] + "-" + el["profile"] for el in self.data["elements"] if el["geometryType"] == "line"]) + ) + ifcMaterialProfileSets = [None for _ in range(len(mpSets))] + for i, mpSet in enumerate(mpSets): + materialIndex = [mat["referenceName"] for mat in self.data["db"]["materials"]].index(mpSet.split("-")[0]) + profileIndex = [prof["referenceName"] for prof in self.data["db"]["profiles"]].index(mpSet.split("-")[1]) + material = ifcMaterials[materialIndex] + profile = ifcProfiles[profileIndex] + matProf = self.f.createIfcMaterialProfile( + self.data["db"]["materials"][materialIndex]["name"] + + " | " + + self.data["db"]["profiles"][profileIndex]["profileName"], + None, + material, + profile, + ) + ifcMaterialProfileSets[i] = self.f.createIfcMaterialProfileSet(None, None, (matProf,)) + + # create structural elements + ifcElements = [None for _ in range(len(self.data["elements"]))] + for i, el in enumerate(self.data["elements"]): + # geometry - product definition shape + prodDefShape = self.create_geometry(el) + + if el["geometryType"] == "line": + # z axis TODO: group by elements + localZAxis = self.f.createIfcDirection(tuple(el["orientation"][2])) + # element + ifcElements[i] = self.f.createIfcStructuralCurveMember( + self.guid(), + ownerHistory, + el["name"], + None, + None, + localPlacement, + prodDefShape, + el["predefinedType"], + localZAxis, + ) + + if el["geometryType"] == "surface": + ifcElements[i] = self.f.createIfcStructuralSurfaceMember( + self.guid(), + ownerHistory, + el["name"], + None, + None, + localPlacement, + prodDefShape, + el["predefinedType"], + el["thickness"], + ) + + # create structural point connections + ifcConnections = [None for _ in range(len(self.data["connections"]))] + for i, conn in enumerate(self.data["connections"]): + # geometry - product definition shape + prodDefShape = self.create_geometry(conn) + + # boundary conditions + if conn["appliedCondition"]: + bc = self.create_applied_conditions(conn["appliedCondition"], conn["geometryType"]) + if conn["geometryType"] == "point": + appliedCondition = self.f.createIfcBoundaryNodeCondition( + None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] + ) + if conn["geometryType"] == "line": + appliedCondition = self.f.createIfcBoundaryEdgeCondition( + None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] + ) + if conn["geometryType"] == "surface": + appliedCondition = self.f.createIfcBoundaryFaceCondition(None, bc["dx"], bc["dy"], bc["dz"]) + else: + appliedCondition = None + + if conn["geometryType"] == "point": + # local axes + localAxes = self.create_orientation(conn["orientation"]) + # connection + ifcConnections[i] = self.f.createIfcStructuralPointConnection( + self.guid(), + ownerHistory, + conn["name"], + None, + None, + localPlacement, + prodDefShape, + appliedCondition, + localAxes, + ) + + if conn["geometryType"] == "line": + # z axis TODO: group by elements + localZAxis = self.f.createIfcDirection(tuple(conn["orientation"][2])) + # connection + ifcConnections[i] = self.f.createIfcStructuralCurveConnection( + self.guid(), + ownerHistory, + conn["name"], + None, + None, + localPlacement, + prodDefShape, + appliedCondition, + localZAxis, + ) + + if conn["geometryType"] == "surface": + ifcConnections[i] = self.f.createIfcStructuralSurfaceConnection( + self.guid(), ownerHistory, conn["name"], None, None, localPlacement, prodDefShape, appliedCondition + ) + + # assign material-profile-sets + for i, mpSet in enumerate(mpSets): + groupOfElements = [] + for j, el in enumerate(self.data["elements"]): + if el["geometryType"] == "line" and el["material"] + "-" + el["profile"] == mpSet: + groupOfElements.append(ifcElements[j]) + + if groupOfElements: + self.f.createIfcRelAssociatesMaterial( + self.guid(), ownerHistory, None, None, tuple(groupOfElements), ifcMaterialProfileSets[i] + ) + + # assign materials + for i, mat in enumerate(self.data["db"]["materials"]): + groupOfElements = [] + for j, el in enumerate(self.data["elements"]): + if el["geometryType"] == "surface" and el["material"] == mat["referenceName"]: + groupOfElements.append(ifcElements[j]) + if groupOfElements: + self.f.createIfcRelAssociatesMaterial( + self.guid(), ownerHistory, None, None, tuple(groupOfElements), ifcMaterials[i] + ) + + # create connections with elements + for i, el in enumerate(self.data["elements"]): + for conn in el["connections"]: + j = [c["referenceName"] for c in self.data["connections"]].index(conn["relatedConnection"]) + geometryType = self.data["connections"][j]["geometryType"] + + if conn["appliedCondition"]: + bc = self.create_applied_conditions(conn["appliedCondition"], geometryType) + if geometryType == "point": + appliedCondition = self.f.createIfcBoundaryNodeCondition( + None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] + ) + if geometryType == "line": + appliedCondition = self.f.createIfcBoundaryEdgeCondition( + None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] + ) + if geometryType == "surface": + appliedCondition = self.f.createIfcBoundaryFaceCondition(None, bc["dx"], bc["dy"], bc["dz"]) + else: + appliedCondition = None + + # local axes + localAxes = self.create_orientation(conn["orientation"]) + + if geometryType == "point": + if not conn["eccentricity"]: + self.f.createIfcRelConnectsStructuralMember( + self.guid(), + ownerHistory, + None, + None, + ifcElements[i], + ifcConnections[j], + appliedCondition, + None, + None, + localAxes, + ) + else: + pointOnElement = self.f.createIfcCartesianPoint(tuple(conn["eccentricity"]["pointOnElement"])) + vector = conn["eccentricity"]["vector"] + connPointEcc = self.f.createIfcConnectionPointEccentricity( + pointOnElement, None, vector[0], vector[1], vector[2] + ) + self.f.createIfcRelConnectsWithEccentricity( + self.guid(), + ownerHistory, + None, + None, + ifcElements[i], + ifcConnections[j], + appliedCondition, + None, + None, + localAxes, + connPointEcc, + ) + + if geometryType in ["line", "surface"]: + self.f.createIfcRelConnectsStructuralMember( + self.guid(), + ownerHistory, + None, + None, + ifcElements[i], + ifcConnections[j], + appliedCondition, + None, + None, + localAxes, + ) + + # assign elements and connections to group + self.f.createIfcRelAssignsToGroup( + self.guid(), ownerHistory, None, None, tuple(ifcElements + ifcConnections), None, model + ) + + # finalize ifc file + self.f.write(self.outputFilename) + + def guid(self): + return ifcopenshell.guid.new() + + def create_header(self): + self.f.wrapped_data.header.file_name.name = os.path.basename(self.outputFilename) + + def create_global_axes(self): + self.xAxis = self.f.createIfcDirection((1.0, 0.0, 0.0)) + self.yAxis = self.f.createIfcDirection((0.0, 1.0, 0.0)) + self.zAxis = self.f.createIfcDirection((0.0, 0.0, 1.0)) + self.origin = self.f.createIfcCartesianPoint((0.0, 0.0, 0.0)) + axes = self.f.createIfcAxis2Placement3D(self.origin, self.zAxis, self.xAxis) + + return axes + + def create_orientation(self, orientation): + xAxis = self.f.createIfcDirection(tuple(orientation[0])) + zAxis = self.f.createIfcDirection(tuple(orientation[2])) + axes = self.f.createIfcAxis2Placement3D(self.origin, zAxis, xAxis) + + return axes + + def create_owner_history(self): + actor = self.f.createIfcActorRole("ENGINEER", None, None) + person = self.f.createIfcPerson("Christovasilis", None, "Ioannis", None, None, None, (actor,)) + organization = self.f.createIfcOrganization( + None, + "IfcOpenShell", + "IfcOpenShell, an open source (LGPL) software library that helps users and software developers to work with the IFC file format.", + ) + p_o = self.f.createIfcPersonAndOrganization(person, organization) + application = self.f.createIfcApplication(organization, "v0.0.x", "IFC2CA", "IFC2CA") + timestamp = int(datetime.now().timestamp()) + ownerHistory = self.f.createIfcOwnerHistory(p_o, application, "READWRITE", None, None, None, None, timestamp) + + return ownerHistory + + def create_reference_subrep(self, globalAxes): + modelRep = self.f.createIfcGeometricRepresentationContext(None, "Model", 3, 1.0e-05, globalAxes, None) + bodySubRep = self.f.createIfcGeometricRepresentationSubContext( + "Body", "Model", None, None, None, None, modelRep, None, "MODEL_VIEW", None + ) + refSubRep = self.f.createIfcGeometricRepresentationSubContext( + "Reference", "Model", None, None, None, None, modelRep, None, "GRAPH_VIEW", None + ) + + return {"model": modelRep, "body": bodySubRep, "reference": refSubRep} + + def create_material(self, material): + ifcMaterial = self.f.createIfcMaterial(material["name"], None, material["category"]) + + mechProps = [] + if "youngModulus" in material["mechProps"]: + youngModulus = self.f.createIfcPropertySingleValue( + "YoungModulus", None, self.f.createIfcModulusOfElasticityMeasure(material["mechProps"]["youngModulus"]) + ) + mechProps.append(youngModulus) + if "shearModulus" in material["mechProps"]: + shearModulus = self.f.createIfcPropertySingleValue( + "ShearModulus", None, self.f.createIfcModulusOfElasticityMeasure(material["mechProps"]["shearModulus"]) + ) + mechProps.append(shearModulus) + if "poissonRatio" in material["mechProps"]: + poissonRatio = self.f.createIfcPropertySingleValue( + "PoissonRatio", None, self.f.createIfcPositiveRatioMeasure(material["mechProps"]["poissonRatio"]) + ) + mechProps.append(poissonRatio) + if mechProps: + self.f.createIfcMaterialProperties( + "Pset_MaterialMechanical", material["name"], tuple(mechProps), ifcMaterial + ) + + commonProps = [] + if "massDensity" in material["commonProps"]: + massDensity = self.f.createIfcPropertySingleValue( + "MassDensity", None, self.f.createIfcMassDensityMeasure(material["commonProps"]["massDensity"]) + ) + commonProps.append(massDensity) + if commonProps: + self.f.createIfcMaterialProperties("Pset_MaterialCommon", material["name"], tuple(commonProps), ifcMaterial) + + return ifcMaterial + + def create_profile(self, profile): + if profile["profileShape"] == "rectangular": + ifcProfile = self.f.createIfcRectangleProfileDef( + profile["profileType"], profile["profileName"], None, profile["xDim"], profile["yDim"] + ) + + if profile["profileShape"] == "iSymmetrical": + ifcProfile = self.f.createIfcIShapeProfileDef( + profile["profileType"], + profile["profileName"], + None, + profile["commonProps"]["overallWidth"], + profile["commonProps"]["overallDepth"], + profile["commonProps"]["webThickness"], + profile["commonProps"]["flangeThickness"], + profile["commonProps"]["filletRadius"], + ) + + mechProps = [] + if "massPerLength" in profile["mechProps"]: + massPerLength = self.f.createIfcPropertySingleValue( + "MassPerLength", None, self.f.createIfcMassPerLengthMeasure(profile["mechProps"]["massPerLength"]) + ) + mechProps.append(massPerLength) + if "crossSectionArea" in profile["mechProps"]: + crossSectionArea = self.f.createIfcPropertySingleValue( + "CrossSectionArea", None, self.f.createIfcAreaMeasure(profile["mechProps"]["crossSectionArea"]) + ) + mechProps.append(crossSectionArea) + if "momentOfInertiaY" in profile["mechProps"]: + momentOfInertiaY = self.f.createIfcPropertySingleValue( + "MomentOfInertiaY", + None, + self.f.createIfcMomentOfInertiaMeasure(profile["mechProps"]["momentOfInertiaY"]), + ) + mechProps.append(momentOfInertiaY) + if "momentOfInertiaZ" in profile["mechProps"]: + momentOfInertiaZ = self.f.createIfcPropertySingleValue( + "MomentOfInertiaZ", + None, + self.f.createIfcMomentOfInertiaMeasure(profile["mechProps"]["momentOfInertiaZ"]), + ) + mechProps.append(momentOfInertiaZ) + if "torsionalConstantX" in profile["mechProps"]: + torsionalConstantX = self.f.createIfcPropertySingleValue( + "TorsionalConstantX", + None, + self.f.createIfcMomentOfInertiaMeasure(profile["mechProps"]["torsionalConstantX"]), + ) + mechProps.append(torsionalConstantX) + if mechProps: + self.f.createIfcProfileProperties( + "Pset_ProfileMechanical", profile["profileName"], tuple(mechProps), ifcProfile + ) + + return ifcProfile + + def create_geometry(self, object): + if object["geometryType"] == "point": + point = self.f.createIfcCartesianPoint(tuple(object["geometry"])) + vertex = self.f.createIfcVertexPoint(point) + vertexTopologyRep = self.f.createIfcTopologyRepresentation( + self.reps["reference"], "Reference", "Vertex", (vertex,) + ) + vertexProdDefShape = self.f.createIfcProductDefinitionShape(None, None, (vertexTopologyRep,)) + + return vertexProdDefShape + + if object["geometryType"] == "line": + startPoint = self.f.createIfcCartesianPoint(tuple(object["geometry"][0])) + startVertex = self.f.createIfcVertexPoint(startPoint) + endPoint = self.f.createIfcCartesianPoint(tuple(object["geometry"][1])) + endVertex = self.f.createIfcVertexPoint(endPoint) + edge = self.f.createIfcEdge(startVertex, endVertex) + edgeTopologyRep = self.f.createIfcTopologyRepresentation( + self.reps["reference"], "Reference", "Edge", (edge,) + ) + edgeProdDefShape = self.f.createIfcProductDefinitionShape(None, None, (edgeTopologyRep,)) + + return edgeProdDefShape + + if object["geometryType"] == "surface": + verts = [None for _ in range(len(object["geometry"]))] + for i, p in enumerate(object["geometry"]): + point = self.f.createIfcCartesianPoint(tuple(p)) + verts[i] = self.f.createIfcVertexPoint(point) + + orientedEdges = [None for _ in range(len(object["geometry"]))] + for i, v in enumerate(verts): + v2Index = (i + 1) if i < len(verts) - 1 else 0 + edge = self.f.createIfcEdge(v, verts[v2Index]) + orientedEdges[i] = self.f.createIfcOrientedEdge(None, None, edge, True) + + edgeLoop = self.f.createIfcEdgeLoop(tuple(orientedEdges)) + localAxes = self.create_orientation(object["orientation"]) + plane = self.f.createIfcPlane(localAxes) + faceBound = self.f.createIfcFaceBound(edgeLoop, True) + face = self.f.createIfcFaceSurface((faceBound,), plane, True) + faceTopologyRep = self.f.createIfcTopologyRepresentation( + self.reps["reference"], "Reference", "Face", (face,) + ) + faceProdDefShape = self.f.createIfcProductDefinitionShape(None, None, (faceTopologyRep,)) + + return faceProdDefShape + + def create_applied_conditions(self, bc, geometryType): + for dof in ["dx", "dy", "dz"]: + if isinstance(bc[dof], bool): + bc[dof] = self.f.createIfcBoolean(bc[dof]) + else: + if geometryType == "point": + bc[dof] = self.f.createIfcLinearStiffnessMeasure(bc[dof]) + if geometryType == "line": + bc[dof] = self.f.createIfcModulusOfLinearSubgradeReactionMeasure(bc[dof]) + if geometryType == "surface": + bc[dof] = self.f.createIfcModulusOfSubgradeReactionMeasure(bc[dof]) + + for dof in ["drx", "dry", "drz"]: + if isinstance(bc[dof], bool): + bc[dof] = self.f.createIfcBoolean(bc[dof]) + else: + if geometryType == "point": + bc[dof] = self.f.createIfcRotationalStiffnessMeasure(bc[dof]) + if geometryType == "line": + bc[dof] = self.f.createIfcModulusOfRotationalSubgradeReactionMeasure(bc[dof]) + + return bc + + +if __name__ == "__main__": + inputFilename = "grid_of_beams.json" + outputFilename = "grid_of_beams.ifc" + + ca2ifc = CA2IFC(inputFilename, outputFilename) + ca2ifc.convert() diff --git a/src/ifc2ca/scriptCodeAsterBonded.py b/src/ifc2ca/_deprecated/scriptCodeAsterBonded.py similarity index 100% rename from src/ifc2ca/scriptCodeAsterBonded.py rename to src/ifc2ca/_deprecated/scriptCodeAsterBonded.py diff --git a/src/ifc2ca/scriptSalomeBonded.py b/src/ifc2ca/_deprecated/scriptSalomeBonded.py similarity index 100% rename from src/ifc2ca/scriptSalomeBonded.py rename to src/ifc2ca/_deprecated/scriptSalomeBonded.py diff --git a/src/ifc2ca/ca2ifc.py b/src/ifc2ca/ca2ifc.py index 7847a88e80..151edce378 100644 --- a/src/ifc2ca/ca2ifc.py +++ b/src/ifc2ca/ca2ifc.py @@ -1,6 +1,6 @@ # Ifc2CA - IFC Code_Aster utility -# Copyright (C) 2020, 2021 Ioannis P. Christovasilis +# Copyright (C) 2020, 2021, 2023, 2024 Ioannis P. Christovasilis # # This file is part of Ifc2CA. # @@ -17,509 +17,435 @@ # You should have received a copy of the GNU Lesser General Public License # along with Ifc2CA. If not, see . -import json -import ifcopenshell -import os -from datetime import datetime +import itertools + +import ifcopenshell as ios +import meshio +import numpy as np + +flatten = itertools.chain.from_iterable -class CA2IFC: - def __init__(self, inputFilename, outputFilename): - self.inputFilename = inputFilename - self.outputFilename = outputFilename - self.data = None - self.f = None - self.reps = {} - self.origin = None - self.xAxis = None - self.yAxis = None - self.zAxis = None +def get_element_data(model, name, element): + if element["geometry_type"] == "Edge": + for i, cell_block in enumerate(model.cells): + if cell_block.type == "line": + cell_tags = model.cell_data["cell_tags"][i] + break + rows = [] + for i_row, i in enumerate(cell_tags): + if i == 0: + continue + tags = model.cell_tags[i] + for tag in tags: + if tag == name: + # print(i_row, i) + rows.append(i_row) + break - def convert(self): - # load json file - with open(self.inputFilename) as dataFile: - self.data = json.load(dataFile) + points = list(set(flatten([cell_block.data[c] for c in rows]))) + points.sort(key=lambda p: np.linalg.norm(model.points[p] - np.array(element["origin"]))) + coords = [np.round(model.points[p], 4).tolist() for p in points] + local_coords = [ + [float(round(np.linalg.norm(model.points[p] - np.array(element["origin"])), 4))] for p in points + ] - # initiate ifc file - self.f = ifcopenshell.file() + return { + "name": name, + "points": points, + "coords": coords, + "local_coords": local_coords, + } - # create header - self.create_header() + elif element["geometry_type"] == "Face": + triangle_cell_tags = None + quad_cell_tags = None + for i, cell_block in enumerate(model.cells): + if cell_block.type == "triangle": + triangle_cell_tags = model.cell_data["cell_tags"][i] + break - # create global axes - globalAxes = self.create_global_axes() - localPlacement = self.f.createIfcLocalPlacement(None, globalAxes) - - # TODO: create units - lengthUnit = self.f.createIfcSIUnit(None, "LENGTHUNIT", None, "METRE") - unitAssignment = self.f.createIfcUnitAssignment((lengthUnit,)) - - # create owner history - ownerHistory = self.create_owner_history() - - # create representations and subrepresentations - self.reps = self.create_reference_subrep(globalAxes) - - # create project and model - project = self.f.createIfcProject( - self.guid(), ownerHistory, "A Project", None, None, None, None, (self.reps["model"],), unitAssignment - ) - model = self.f.createIfcStructuralAnalysisModel( - self.guid(), - ownerHistory, - self.data["name"], - None, - None, - "NOTDEFINED", - globalAxes, - None, - None, - localPlacement, - ) - self.f.createIfcRelDeclares(self.guid(), ownerHistory, None, None, project, (model,)) - - # create materials - ifcMaterials = [None for _ in range(len(self.data["db"]["materials"]))] - for i, material in enumerate(self.data["db"]["materials"]): - ifcMaterials[i] = self.create_material(material) - - # create profiles - ifcProfiles = [None for _ in range(len(self.data["db"]["profiles"]))] - for i, profile in enumerate(self.data["db"]["profiles"]): - ifcProfiles[i] = self.create_profile(profile) - - # create material-profile sets - mpSets = list( - set([el["material"] + "-" + el["profile"] for el in self.data["elements"] if el["geometryType"] == "line"]) - ) - ifcMaterialProfileSets = [None for _ in range(len(mpSets))] - for i, mpSet in enumerate(mpSets): - materialIndex = [mat["referenceName"] for mat in self.data["db"]["materials"]].index(mpSet.split("-")[0]) - profileIndex = [prof["referenceName"] for prof in self.data["db"]["profiles"]].index(mpSet.split("-")[1]) - material = ifcMaterials[materialIndex] - profile = ifcProfiles[profileIndex] - matProf = self.f.createIfcMaterialProfile( - self.data["db"]["materials"][materialIndex]["name"] - + " | " - + self.data["db"]["profiles"][profileIndex]["profileName"], - None, - material, - profile, - ) - ifcMaterialProfileSets[i] = self.f.createIfcMaterialProfileSet(None, None, (matProf,)) - - # create structural elements - ifcElements = [None for _ in range(len(self.data["elements"]))] - for i, el in enumerate(self.data["elements"]): - # geometry - product definition shape - prodDefShape = self.create_geometry(el) - - if el["geometryType"] == "line": - # z axis TODO: group by elements - localZAxis = self.f.createIfcDirection(tuple(el["orientation"][2])) - # element - ifcElements[i] = self.f.createIfcStructuralCurveMember( - self.guid(), - ownerHistory, - el["name"], - None, - None, - localPlacement, - prodDefShape, - el["predefinedType"], - localZAxis, - ) - - if el["geometryType"] == "surface": - ifcElements[i] = self.f.createIfcStructuralSurfaceMember( - self.guid(), - ownerHistory, - el["name"], - None, - None, - localPlacement, - prodDefShape, - el["predefinedType"], - el["thickness"], - ) - - # create structural point connections - ifcConnections = [None for _ in range(len(self.data["connections"]))] - for i, conn in enumerate(self.data["connections"]): - # geometry - product definition shape - prodDefShape = self.create_geometry(conn) - - # boundary conditions - if conn["appliedCondition"]: - bc = self.create_applied_conditions(conn["appliedCondition"], conn["geometryType"]) - if conn["geometryType"] == "point": - appliedCondition = self.f.createIfcBoundaryNodeCondition( - None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] - ) - if conn["geometryType"] == "line": - appliedCondition = self.f.createIfcBoundaryEdgeCondition( - None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] - ) - if conn["geometryType"] == "surface": - appliedCondition = self.f.createIfcBoundaryFaceCondition(None, bc["dx"], bc["dy"], bc["dz"]) + if triangle_cell_tags is not None: + rows = [] + for i_row, i in enumerate(triangle_cell_tags): + if i == 0: + continue + tags = model.cell_tags[i] + for tag in tags: + if tag == name: + # print(i_row, i) + rows.append(i_row) + break + if not len(rows): + points = [] else: - appliedCondition = None + points = list(flatten([cell_block.data[c] for c in rows])) - if conn["geometryType"] == "point": - # local axes - localAxes = self.create_orientation(conn["orientation"]) - # connection - ifcConnections[i] = self.f.createIfcStructuralPointConnection( - self.guid(), - ownerHistory, - conn["name"], - None, - None, - localPlacement, - prodDefShape, - appliedCondition, - localAxes, - ) + for i, cell_block in enumerate(model.cells): + if cell_block.type == "quad": + quad_cell_tags = model.cell_data["cell_tags"][i] + break - if conn["geometryType"] == "line": - # z axis TODO: group by elements - localZAxis = self.f.createIfcDirection(tuple(conn["orientation"][2])) - # connection - ifcConnections[i] = self.f.createIfcStructuralCurveConnection( - self.guid(), - ownerHistory, - conn["name"], - None, - None, - localPlacement, - prodDefShape, - appliedCondition, - localZAxis, - ) + if quad_cell_tags is not None: + rows = [] + for i_row, i in enumerate(quad_cell_tags): + if i == 0: + continue + tags = model.cell_tags[i] + for tag in tags: + if tag == name: + # print(i_row, i) + rows.append(i_row) + break + if len(rows): + points.extend(list(flatten([cell_block.data[c] for c in rows]))) - if conn["geometryType"] == "surface": - ifcConnections[i] = self.f.createIfcStructuralSurfaceConnection( - self.guid(), ownerHistory, conn["name"], None, None, localPlacement, prodDefShape, appliedCondition - ) + points = list(set(points)) + points.sort() + coords = [model.points[p].tolist() for p in points] + local_coords = [ + np.round(np.array(element["orientation"]).dot(model.points[p] - np.array(element["origin"])), 4).tolist()[ + :2 + ] + for p in points + ] - # assign material-profile-sets - for i, mpSet in enumerate(mpSets): - groupOfElements = [] - for j, el in enumerate(self.data["elements"]): - if el["geometryType"] == "line" and el["material"] + "-" + el["profile"] == mpSet: - groupOfElements.append(ifcElements[j]) + return { + "name": name, + "points": points, + "coords": coords, + "local_coords": local_coords, + } - if groupOfElements: - self.f.createIfcRelAssociatesMaterial( - self.guid(), ownerHistory, None, None, tuple(groupOfElements), ifcMaterialProfileSets[i] - ) - # assign materials - for i, mat in enumerate(self.data["db"]["materials"]): - groupOfElements = [] - for j, el in enumerate(self.data["elements"]): - if el["geometryType"] == "surface" and el["material"] == mat["referenceName"]: - groupOfElements.append(ifcElements[j]) - if groupOfElements: - self.f.createIfcRelAssociatesMaterial( - self.guid(), ownerHistory, None, None, tuple(groupOfElements), ifcMaterials[i] - ) +def get_element_result_data(model, field_label, name, element, field_type): + points = get_element_data(model, name, element)["points"] + if field_type == "InternalForces": + if element["geometry_type"] == "Edge": + return { + "N": [round(model.point_data[field_label][p][0], 4) for p in points], + "VY": [round(model.point_data[field_label][p][1], 4) for p in points], + "VZ": [round(model.point_data[field_label][p][2], 4) for p in points], + "MT": [round(model.point_data[field_label][p][3], 4) for p in points], + "MFY": [round(model.point_data[field_label][p][4], 4) for p in points], + "MFZ": [round(model.point_data[field_label][p][5], 4) for p in points], + } - # create connections with elements - for i, el in enumerate(self.data["elements"]): - for conn in el["connections"]: - j = [c["referenceName"] for c in self.data["connections"]].index(conn["relatedConnection"]) - geometryType = self.data["connections"][j]["geometryType"] + elif element["geometry_type"] == "Face": + if len(model.point_data[field_label][points[0]]) == 8: + offset = 0 + elif len(model.point_data[field_label][points[0]]) == 14: + offset = 6 + else: + assert ( + False + ), f"Internal force field with {len(model.point_data[field_label][points[0]])} field values for {field_label} and {element['Name']} " - if conn["appliedCondition"]: - bc = self.create_applied_conditions(conn["appliedCondition"], geometryType) - if geometryType == "point": - appliedCondition = self.f.createIfcBoundaryNodeCondition( - None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] - ) - if geometryType == "line": - appliedCondition = self.f.createIfcBoundaryEdgeCondition( - None, bc["dx"], bc["dy"], bc["dz"], bc["drx"], bc["dry"], bc["drz"] - ) - if geometryType == "surface": - appliedCondition = self.f.createIfcBoundaryFaceCondition(None, bc["dx"], bc["dy"], bc["dz"]) - else: - appliedCondition = None + return { + "NXX": [round(model.point_data[field_label][p][offset + 0], 4) for p in points], + "NYY": [round(model.point_data[field_label][p][offset + 1], 4) for p in points], + "NXY": [round(model.point_data[field_label][p][offset + 2], 4) for p in points], + "MXX": [round(model.point_data[field_label][p][offset + 3], 4) for p in points], + "MYY": [round(model.point_data[field_label][p][offset + 4], 4) for p in points], + "MXY": [round(model.point_data[field_label][p][offset + 5], 4) for p in points], + "QX": [round(model.point_data[field_label][p][offset + 6], 4) for p in points], + "QY": [round(model.point_data[field_label][p][offset + 7], 4) for p in points], + } - # local axes - localAxes = self.create_orientation(conn["orientation"]) + if field_type == "Displacements": + return { + "DX": [round(model.point_data[field_label][p][0], 4) for p in points], + "DY": [round(model.point_data[field_label][p][1], 4) for p in points], + "DZ": [round(model.point_data[field_label][p][2], 4) for p in points], + "DRX": [round(model.point_data[field_label][p][3], 4) for p in points], + "DRY": [round(model.point_data[field_label][p][4], 4) for p in points], + "DRZ": [round(model.point_data[field_label][p][5], 4) for p in points], + } - if geometryType == "point": - if not conn["eccentricity"]: - self.f.createIfcRelConnectsStructuralMember( - self.guid(), - ownerHistory, - None, - None, - ifcElements[i], - ifcConnections[j], - appliedCondition, - None, - None, - localAxes, - ) - else: - pointOnElement = self.f.createIfcCartesianPoint(tuple(conn["eccentricity"]["pointOnElement"])) - vector = conn["eccentricity"]["vector"] - connPointEcc = self.f.createIfcConnectionPointEccentricity( - pointOnElement, None, vector[0], vector[1], vector[2] - ) - self.f.createIfcRelConnectsWithEccentricity( - self.guid(), - ownerHistory, - None, - None, - ifcElements[i], - ifcConnections[j], - appliedCondition, - None, - None, - localAxes, - connPointEcc, - ) - if geometryType in ["line", "surface"]: - self.f.createIfcRelConnectsStructuralMember( - self.guid(), - ownerHistory, - None, - None, - ifcElements[i], - ifcConnections[j], - appliedCondition, - None, - None, - localAxes, +def results_to_ifc(ifc_file, ifc_model, rmed_path, global_case, field_types, data): + if not rmed_path.exists(): + print(f"Med file with results not found for case_instant: {global_case}") + return + + result = meshio.read(rmed_path, "med") + if global_case == "LC": + model_cases = data["load_cases"] + elif global_case == "COMB": + model_cases = data["load_combinations"] + for field in field_types: + if field == "InternalForces": + _parsed_data = internal_forces_to_ifc(ifc_file, ifc_model, result, model_cases, data["elements"]) + elif field == "Displacements": + _parsed_data = displacements_to_ifc(ifc_file, ifc_model, result, model_cases, data["elements"]) + + +def internal_forces_to_ifc(ifc_file, ifc_model, result, model_cases, elements): + result_cases = [dict() for _ in model_cases] + field_cases = [f"ELEMENT_FORCE[{i}] - {i + 1}" for i in range(len(result_cases))] + + # Create Result Groups for load case_instance combinations + for iCase, case_instance in enumerate(model_cases): + result_cases[iCase]["case_instance"] = ifc_file.create_entity( + "IfcStructuralResultGroup", + **{ + "GlobalId": ios.guid.new(), + "Name": "Internal Forces for " + case_instance["Name"], + "TheoryType": "FIRST_ORDER_THEORY", + "ResultForLoadGroup": ifc_file.by_id(case_instance["id"]), + "IsLinear": True, + }, + ) + + result_cases[iCase]["assignment"] = ifc_file.create_entity( + "IfcRelAssignsToGroup", + **{ + "GlobalId": ios.guid.new(), + "RelatedObjects": [], + "RelatingGroup": result_cases[iCase]["case_instance"], + }, + ) + + if ifc_model.HasResults: + ifc_model.HasResults += tuple([result["case_instance"] for result in result_cases]) + else: + ifc_model.HasResults = tuple([result["case_instance"] for result in result_cases]) + + data = [] + for _, element in enumerate(elements): + group_name = getGroupName(element["ref_id"]) + name = element["Name"] + info = get_element_data(result, group_name, element) + assert len(info["coords"]) >= 2 + for iCase, field_case in enumerate(field_cases): + forces = get_element_result_data(result, field_case, group_name, element, field_type="InternalForces") + reaction = ifc_file.create_entity( + "IfcStructuralCurveReaction" if element["geometry_type"] == "Edge" else "IfcStructuralSurfaceReaction", + **{ + "GlobalId": ios.guid.new(), + "Name": "Internal Forces for " + model_cases[iCase]["Name"] + f" on {name}", + # "AppliedLoad": load["ifcLoad"], + "GlobalOrLocal": "LOCAL_COORDS", + "PredefinedType": "DISCRETE", + }, + ) + result_cases[iCase]["assignment"].RelatedObjects += (reaction,) + + ifc_file.create_entity( + "IfcRelConnectsStructuralActivity", + **{ + "GlobalId": ios.guid.new(), + "RelatingElement": ifc_file.by_id(element["id"]), + "RelatedStructuralActivity": reaction, + }, + ) + + reaction.AppliedLoad = ifc_file.create_entity( + "IfcStructuralLoadConfiguration", + **{ + "Name": "Internal Forces for " + model_cases[iCase]["Name"] + f" on {name}", + "Values": [], + "Locations": tuple([tuple(node) for node in info["local_coords"]]), + }, + ) + + if element["geometry_type"] == "Edge": + for iNode, node in enumerate(info["coords"]): + location = f"({node[0]}, {node[1]}, {node[2]})" + distance = info["local_coords"][iNode][0] + + N = forces["N"][iNode] + VY = forces["VY"][iNode] + VZ = forces["VZ"][iNode] + MT = forces["MT"][iNode] + MFY = forces["MFY"][iNode] + MFZ = forces["MFZ"][iNode] + + data.append([name, f"LCC-{iCase + 1} @ {distance}", location, N, VY, VZ, MT, MFY, MFZ]) + + pointValue = ifc_file.create_entity( + "IfcStructuralLoadSingleForce", + **{ + "Name": "Internal Forces for " + model_cases[iCase]["Name"] + f" @ {distance} on {name}", + "ForceX": N, + "ForceY": VY, + "ForceZ": VZ, + "MomentX": MT, + "MomentY": MFY, + "MomentZ": MFZ, + }, ) + reaction.AppliedLoad.Values += (pointValue,) - # assign elements and connections to group - self.f.createIfcRelAssignsToGroup( - self.guid(), ownerHistory, None, None, tuple(ifcElements + ifcConnections), None, model + elif element["geometry_type"] == "Face": + for iNode, node in enumerate(info["coords"]): + location = f"({node[0]}, {node[1]}, {node[2]})" + distance = tuple(info["local_coords"][iNode]) + + NXX = forces["NXX"][iNode] + NYY = forces["NYY"][iNode] + NXY = forces["NXY"][iNode] + MXX = forces["MXX"][iNode] + MYY = forces["MYY"][iNode] + MXY = forces["MXY"][iNode] + + data.append([name, f"LCC-{iCase + 1} @ {distance}", location, NXX, NYY, NXY, MXX, MYY, MXY]) + + pointValue = ifc_file.create_entity( + "IfcStructuralLoadSingleForce", + **{ + "Name": "Internal Forces for " + model_cases[iCase]["Name"] + f" @ {distance} on {name}", + "ForceX": NXX, + "ForceY": NYY, + "ForceZ": NXY, + "MomentX": MXX, + "MomentY": MYY, + "MomentZ": MXY, + }, + ) + reaction.AppliedLoad.Values += (pointValue,) + + return data + + +def displacements_to_ifc(ifc_file, ifc_model, result, model_cases, elements): + result_cases = [dict() for _ in model_cases] + field_cases = [f"MODEL_DISP[{i}] - {i + 1}" for i in range(len(result_cases))] + + # Create Result Groups for load case_instance combinations + for iCase, case_instance in enumerate(model_cases): + result_cases[iCase]["case_instance"] = ifc_file.create_entity( + "IfcStructuralResultGroup", + **{ + "GlobalId": ios.guid.new(), + "Name": "Global Displacements for " + case_instance["Name"], + "TheoryType": "FIRST_ORDER_THEORY", + "ResultForLoadGroup": ifc_file.by_id(case_instance["id"]), + "IsLinear": True, + }, ) - # finalize ifc file - self.f.write(self.outputFilename) - - def guid(self): - return ifcopenshell.guid.new() - - def create_header(self): - self.f.wrapped_data.header.file_name.name = os.path.basename(self.outputFilename) - - def create_global_axes(self): - self.xAxis = self.f.createIfcDirection((1.0, 0.0, 0.0)) - self.yAxis = self.f.createIfcDirection((0.0, 1.0, 0.0)) - self.zAxis = self.f.createIfcDirection((0.0, 0.0, 1.0)) - self.origin = self.f.createIfcCartesianPoint((0.0, 0.0, 0.0)) - axes = self.f.createIfcAxis2Placement3D(self.origin, self.zAxis, self.xAxis) - - return axes - - def create_orientation(self, orientation): - xAxis = self.f.createIfcDirection(tuple(orientation[0])) - zAxis = self.f.createIfcDirection(tuple(orientation[2])) - axes = self.f.createIfcAxis2Placement3D(self.origin, zAxis, xAxis) - - return axes - - def create_owner_history(self): - actor = self.f.createIfcActorRole("ENGINEER", None, None) - person = self.f.createIfcPerson("Christovasilis", None, "Ioannis", None, None, None, (actor,)) - organization = self.f.createIfcOrganization( - None, - "IfcOpenShell", - "IfcOpenShell, an open source (LGPL) software library that helps users and software developers to work with the IFC file format.", - ) - p_o = self.f.createIfcPersonAndOrganization(person, organization) - application = self.f.createIfcApplication(organization, "v0.0.x", "IFC2CA", "IFC2CA") - timestamp = int(datetime.now().timestamp()) - ownerHistory = self.f.createIfcOwnerHistory(p_o, application, "READWRITE", None, None, None, None, timestamp) - - return ownerHistory - - def create_reference_subrep(self, globalAxes): - modelRep = self.f.createIfcGeometricRepresentationContext(None, "Model", 3, 1.0e-05, globalAxes, None) - bodySubRep = self.f.createIfcGeometricRepresentationSubContext( - "Body", "Model", None, None, None, None, modelRep, None, "MODEL_VIEW", None - ) - refSubRep = self.f.createIfcGeometricRepresentationSubContext( - "Reference", "Model", None, None, None, None, modelRep, None, "GRAPH_VIEW", None + result_cases[iCase]["assignment"] = ifc_file.create_entity( + "IfcRelAssignsToGroup", + **{ + "GlobalId": ios.guid.new(), + "RelatedObjects": [], + "RelatingGroup": result_cases[iCase]["case_instance"], + }, ) - return {"model": modelRep, "body": bodySubRep, "reference": refSubRep} + if ifc_model.HasResults: + ifc_model.HasResults += tuple([result["case_instance"] for result in result_cases]) + else: + ifc_model.HasResults = tuple([result["case_instance"] for result in result_cases]) - def create_material(self, material): - ifcMaterial = self.f.createIfcMaterial(material["name"], None, material["category"]) + data = [] + for _, element in enumerate(elements): + group_name = getGroupName(element["ref_id"]) + name = element["Name"] + info = get_element_data(result, group_name, element) + assert len(info["coords"]) >= 2 + for iCase, case_instance in enumerate(field_cases): + displacements = get_element_result_data( + result, case_instance, group_name, element, field_type="Displacements" + ) + reaction = ifc_file.create_entity( + "IfcStructuralCurveReaction" if element["geometry_type"] == "Edge" else "IfcStructuralSurfaceReaction", + **{ + "GlobalId": ios.guid.new(), + "Name": "Global Displacements for " + model_cases[iCase]["Name"] + f" on {name}", + # "AppliedLoad": load["ifcLoad"], + "GlobalOrLocal": "LOCAL_COORDS", + "PredefinedType": "DISCRETE", + }, + ) + result_cases[iCase]["assignment"].RelatedObjects += (reaction,) - mechProps = [] - if "youngModulus" in material["mechProps"]: - youngModulus = self.f.createIfcPropertySingleValue( - "YoungModulus", None, self.f.createIfcModulusOfElasticityMeasure(material["mechProps"]["youngModulus"]) - ) - mechProps.append(youngModulus) - if "shearModulus" in material["mechProps"]: - shearModulus = self.f.createIfcPropertySingleValue( - "ShearModulus", None, self.f.createIfcModulusOfElasticityMeasure(material["mechProps"]["shearModulus"]) - ) - mechProps.append(shearModulus) - if "poissonRatio" in material["mechProps"]: - poissonRatio = self.f.createIfcPropertySingleValue( - "PoissonRatio", None, self.f.createIfcPositiveRatioMeasure(material["mechProps"]["poissonRatio"]) - ) - mechProps.append(poissonRatio) - if mechProps: - self.f.createIfcMaterialProperties( - "Pset_MaterialMechanical", material["name"], tuple(mechProps), ifcMaterial + ifc_file.create_entity( + "IfcRelConnectsStructuralActivity", + **{ + "GlobalId": ios.guid.new(), + "RelatingElement": ifc_file.by_id(element["id"]), + "RelatedStructuralActivity": reaction, + }, ) - commonProps = [] - if "massDensity" in material["commonProps"]: - massDensity = self.f.createIfcPropertySingleValue( - "MassDensity", None, self.f.createIfcMassDensityMeasure(material["commonProps"]["massDensity"]) - ) - commonProps.append(massDensity) - if commonProps: - self.f.createIfcMaterialProperties("Pset_MaterialCommon", material["name"], tuple(commonProps), ifcMaterial) - - return ifcMaterial - - def create_profile(self, profile): - if profile["profileShape"] == "rectangular": - ifcProfile = self.f.createIfcRectangleProfileDef( - profile["profileType"], profile["profileName"], None, profile["xDim"], profile["yDim"] + reaction.AppliedLoad = ifc_file.create_entity( + "IfcStructuralLoadConfiguration", + **{ + "Name": "Global Displacements for " + model_cases[iCase]["Name"] + f" on {name}", + "Values": [], + "Locations": tuple([tuple(node) for node in info["local_coords"]]), + }, ) - if profile["profileShape"] == "iSymmetrical": - ifcProfile = self.f.createIfcIShapeProfileDef( - profile["profileType"], - profile["profileName"], - None, - profile["commonProps"]["overallWidth"], - profile["commonProps"]["overallDepth"], - profile["commonProps"]["webThickness"], - profile["commonProps"]["flangeThickness"], - profile["commonProps"]["filletRadius"], - ) + if element["geometry_type"] == "Edge": + for iNode, node in enumerate(info["coords"]): + location = f"({node[0]}, {node[1]}, {node[2]})" + distance = info["local_coords"][iNode][0] - mechProps = [] - if "massPerLength" in profile["mechProps"]: - massPerLength = self.f.createIfcPropertySingleValue( - "MassPerLength", None, self.f.createIfcMassPerLengthMeasure(profile["mechProps"]["massPerLength"]) - ) - mechProps.append(massPerLength) - if "crossSectionArea" in profile["mechProps"]: - crossSectionArea = self.f.createIfcPropertySingleValue( - "CrossSectionArea", None, self.f.createIfcAreaMeasure(profile["mechProps"]["crossSectionArea"]) - ) - mechProps.append(crossSectionArea) - if "momentOfInertiaY" in profile["mechProps"]: - momentOfInertiaY = self.f.createIfcPropertySingleValue( - "MomentOfInertiaY", - None, - self.f.createIfcMomentOfInertiaMeasure(profile["mechProps"]["momentOfInertiaY"]), - ) - mechProps.append(momentOfInertiaY) - if "momentOfInertiaZ" in profile["mechProps"]: - momentOfInertiaZ = self.f.createIfcPropertySingleValue( - "MomentOfInertiaZ", - None, - self.f.createIfcMomentOfInertiaMeasure(profile["mechProps"]["momentOfInertiaZ"]), - ) - mechProps.append(momentOfInertiaZ) - if "torsionalConstantX" in profile["mechProps"]: - torsionalConstantX = self.f.createIfcPropertySingleValue( - "TorsionalConstantX", - None, - self.f.createIfcMomentOfInertiaMeasure(profile["mechProps"]["torsionalConstantX"]), - ) - mechProps.append(torsionalConstantX) - if mechProps: - self.f.createIfcProfileProperties( - "Pset_ProfileMechanical", profile["profileName"], tuple(mechProps), ifcProfile - ) + DX = displacements["DX"][iNode] + DY = displacements["DY"][iNode] + DZ = displacements["DZ"][iNode] + DRX = displacements["DRX"][iNode] + DRY = displacements["DRY"][iNode] + DRZ = displacements["DRZ"][iNode] - return ifcProfile + data.append([name, f"LCC-{iCase + 1} @ {distance}", location, DX, DY, DZ, DRX, DRY, DRZ]) - def create_geometry(self, object): - if object["geometryType"] == "point": - point = self.f.createIfcCartesianPoint(tuple(object["geometry"])) - vertex = self.f.createIfcVertexPoint(point) - vertexTopologyRep = self.f.createIfcTopologyRepresentation( - self.reps["reference"], "Reference", "Vertex", (vertex,) - ) - vertexProdDefShape = self.f.createIfcProductDefinitionShape(None, None, (vertexTopologyRep,)) + pointValue = ifc_file.create_entity( + "IfcStructuralLoadSingleDisplacement", + **{ + "Name": "Global Displacements for " + + model_cases[iCase]["Name"] + + f" @ {distance} on {name}", + "DisplacementX": DX, + "DisplacementY": DY, + "DisplacementZ": DZ, + "RotationalDisplacementRX": DRX, + "RotationalDisplacementRY": DRY, + "RotationalDisplacementRZ": DRZ, + }, + ) + reaction.AppliedLoad.Values += (pointValue,) - return vertexProdDefShape + elif element["geometry_type"] == "Face": + for iNode, node in enumerate(info["coords"]): + location = f"({node[0]}, {node[1]}, {node[2]})" + distance = tuple(info["local_coords"][iNode]) - if object["geometryType"] == "line": - startPoint = self.f.createIfcCartesianPoint(tuple(object["geometry"][0])) - startVertex = self.f.createIfcVertexPoint(startPoint) - endPoint = self.f.createIfcCartesianPoint(tuple(object["geometry"][1])) - endVertex = self.f.createIfcVertexPoint(endPoint) - edge = self.f.createIfcEdge(startVertex, endVertex) - edgeTopologyRep = self.f.createIfcTopologyRepresentation( - self.reps["reference"], "Reference", "Edge", (edge,) - ) - edgeProdDefShape = self.f.createIfcProductDefinitionShape(None, None, (edgeTopologyRep,)) + DX = displacements["DX"][iNode] + DY = displacements["DY"][iNode] + DZ = displacements["DZ"][iNode] + DRX = displacements["DRX"][iNode] + DRY = displacements["DRY"][iNode] + DRZ = displacements["DRZ"][iNode] - return edgeProdDefShape + data.append([name, f"LCC-{iCase + 1} @ {distance}", location, DX, DY, DZ, DRX, DRY, DRZ]) - if object["geometryType"] == "surface": - verts = [None for _ in range(len(object["geometry"]))] - for i, p in enumerate(object["geometry"]): - point = self.f.createIfcCartesianPoint(tuple(p)) - verts[i] = self.f.createIfcVertexPoint(point) + pointValue = ifc_file.create_entity( + "IfcStructuralLoadSingleDisplacement", + **{ + "Name": "Global Displacements for " + + model_cases[iCase]["Name"] + + f" @ {distance} on {name}", + "DisplacementX": DX, + "DisplacementY": DY, + "DisplacementZ": DZ, + "RotationalDisplacementRX": DRX, + "RotationalDisplacementRY": DRY, + "RotationalDisplacementRZ": DRZ, + }, + ) + reaction.AppliedLoad.Values += (pointValue,) - orientedEdges = [None for _ in range(len(object["geometry"]))] - for i, v in enumerate(verts): - v2Index = (i + 1) if i < len(verts) - 1 else 0 - edge = self.f.createIfcEdge(v, verts[v2Index]) - orientedEdges[i] = self.f.createIfcOrientedEdge(None, None, edge, True) - - edgeLoop = self.f.createIfcEdgeLoop(tuple(orientedEdges)) - localAxes = self.create_orientation(object["orientation"]) - plane = self.f.createIfcPlane(localAxes) - faceBound = self.f.createIfcFaceBound(edgeLoop, True) - face = self.f.createIfcFaceSurface((faceBound,), plane, True) - faceTopologyRep = self.f.createIfcTopologyRepresentation( - self.reps["reference"], "Reference", "Face", (face,) - ) - faceProdDefShape = self.f.createIfcProductDefinitionShape(None, None, (faceTopologyRep,)) - - return faceProdDefShape - - def create_applied_conditions(self, bc, geometryType): - for dof in ["dx", "dy", "dz"]: - if isinstance(bc[dof], bool): - bc[dof] = self.f.createIfcBoolean(bc[dof]) - else: - if geometryType == "point": - bc[dof] = self.f.createIfcLinearStiffnessMeasure(bc[dof]) - if geometryType == "line": - bc[dof] = self.f.createIfcModulusOfLinearSubgradeReactionMeasure(bc[dof]) - if geometryType == "surface": - bc[dof] = self.f.createIfcModulusOfSubgradeReactionMeasure(bc[dof]) - - for dof in ["drx", "dry", "drz"]: - if isinstance(bc[dof], bool): - bc[dof] = self.f.createIfcBoolean(bc[dof]) - else: - if geometryType == "point": - bc[dof] = self.f.createIfcRotationalStiffnessMeasure(bc[dof]) - if geometryType == "line": - bc[dof] = self.f.createIfcModulusOfRotationalSubgradeReactionMeasure(bc[dof]) - - return bc + return data -if __name__ == "__main__": - inputFilename = "grid_of_beams.json" - outputFilename = "grid_of_beams.ifc" - - ca2ifc = CA2IFC(inputFilename, outputFilename) - ca2ifc.convert() +def getGroupName(name): + if "|" in name: + info = name.split("|") + sortName = "".join(c for c in info[0] if c.isupper()) + return f"{sortName[2:]}_{info[1]}" + else: + return name diff --git a/src/ifc2ca/ifc2ca.py b/src/ifc2ca/ifc2ca.py index 9231329ee3..2dd4b04263 100644 --- a/src/ifc2ca/ifc2ca.py +++ b/src/ifc2ca/ifc2ca.py @@ -1,5 +1,5 @@ # Ifc2CA - IFC Code_Aster utility -# Copyright (C) 2020, 2021 Ioannis P. Christovasilis +# Copyright (C) 2020, 2021, 2023, 2024 Ioannis P. Christovasilis # # This file is part of Ifc2CA. # @@ -16,618 +16,944 @@ # You should have received a copy of the GNU Lesser General Public License # along with Ifc2CA. If not, see . -from __future__ import division -from __future__ import print_function import json -import ifcopenshell -import numpy as np +import os +import subprocess +from copy import deepcopy from pathlib import Path +from typing import Dict, List + +import ifcopenshell as ios +# import ifcopenshell.geom +import ifcopenshell.util.element +import ifcopenshell.util.placement +import ifcopenshell.util.representation +# import ifcopenshell.util.shape +import numpy as np +from jinja2 import Environment, FileSystemLoader + +from . import ca2ifc +from .scriptCodeAster import CommandFileConstructor -class IFC2CA: - def __init__(self, filename): - self.filename = filename - self.file = None - self.result = {} - self.warnings = [] +class Ifc2CA: + folder_path = None + salome_path = None + model_keys = ["id", "type", "GlobalId", "Name", "LoadedBy", "HasResults"] + member_keys = ["id", "type", "GlobalId", "Name", "Thickness"] + connection_keys = ["id", "type", "GlobalId", "Name"] + material_keys = ["id", "type", "Name", "Category"] + material_property_keys = ["MassDensity", "YoungModulus", "PoissonRatio", "ShearModulus"] + profile_property_keys = ["CrossSectionArea", "MomentOfInertiaY", "MomentOfInertiaZ", "TorsionalConstantX"] + fixed_conditions = { + "Vertex": dict(zip(["dx", "dy", "dz", "drx", "dry", "drz"], [True, True, True, True, True, True])), + "Edge": dict(zip(["dx", "dy", "dz", "drx", "dry", "drz"], [True, True, True, True, True, True])), + "Face": dict(zip(["dx", "dy", "dz"], [True, True, True])), + } + + def __init__(self, path: os.PathLike | str): + self.path = Path(path) self.tol = 1e-06 + self.file = ios.open(self.path) + self.folder_path = self.path.parent / f"{self.path.stem}_ifc2ca" + self.env = Environment(loader=FileSystemLoader(Path(__file__).parent / "templates")) - def convert(self): - self.file = ifcopenshell.open(self.filename) - for model in self.file.by_type("IfcStructuralAnalysisModel"): - elements = self.get_structural_items(model, item_type="IfcStructuralMember") - connections = self.get_structural_items( - model, item_type="IfcStructuralConnection" - ) + # expose GET api functions + ## functions related to models + def get_models(self): + return self.file.by_type("IfcStructuralAnalysisModel") - materialdb = [] - materials = list(dict.fromkeys([e["material"] for e in elements])) - for mat in [mat for mat in materials if mat]: - id = int(mat.split("|")[1]) - material = self.get_material_properties(self.file.by_id(id)) - material["relatedElements"] = [ - e["referenceName"] - for e in elements - if "material" in e and e["material"] == mat - ] - materialdb.append(material) + def get_context(self): + return ifcopenshell.util.representation.get_context(self.file, "Model", "Reference", "GRAPH_VIEW") - profiledb = [] - profiles = list( - dict.fromkeys([e["profile"] for e in elements if "profile" in e]) - ) - for prof in [prof for prof in profiles if prof]: - id = int(prof.split("|")[1]) - profile = self.get_profile_properties(self.file.by_id(id)) - profile["relatedElements"] = [ - e["referenceName"] - for e in elements - if "profile" in e and e["profile"] == prof - ] - profiledb.append(profile) + ## functions related to members and connections + def get_items(self, model: ios.entity_instance | None = None): + if model is not None: + return ifcopenshell.util.element.get_grouped_by(model) + return self.file.by_type("IfcStructuralItem") - self.result = { - "referenceName": model.is_a() + "|" + str(model.id()), - "name": model.Name, - "id": model.GlobalId, - "elements": elements, - "connections": connections, - "db": {"materials": materialdb, "profiles": profiledb}, - "warnings": self.warnings, - } - - print(f"Model {model.Name} converted") - print(f"Number of elements: {len(elements)}") - print(f"Number of connections: {len(connections)}") - print(f"Number of materials: {len(materialdb)}") - print(f"Number of profiles: {len(profiledb)}") - print("") - - break - - def get_structural_items(self, model, item_type="IfcStructuralItem"): - items = [] - for group in model.IsGroupedBy: - for item in group.RelatedObjects: - if not item.is_a(item_type): - continue - data = self.get_item_data(item) - if data: - items.append(data) - return items - - def get_item_data(self, item): - transformation = self.get_transformation(item.ObjectPlacement) - - if item.is_a("IfcStructuralCurveMember"): - representation = self.get_representation(item, "Edge") - material_profile = self.get_material_profile(item) - if not representation: - self.warnings.append( - f"No representation defined for {item.is_a()}|{item.id()}. Member excluded" - ) - return - if not material_profile: - self.warnings.append(f"No material defined for {item.is_a()}|{item.id()}") - self.warnings.append(f"No profile defined for {item.is_a()}|{item.id()}") - materialId = None - profileId = None - else: - material = material_profile.Material - materialId = material.is_a() + "|" + str(material.id()) - profile = material_profile.Profile - profileId = profile.is_a() + "|" + str(profile.id()) - - geometry = self.get_geometry(representation) - orientation = self.get_1D_orientation(geometry, item.Axis) - connections = self.get_connection_data(item.ConnectedBy) - for conn in connections: - if not conn["orientation"]: - conn["orientation"] = orientation - # --> Correct pointOnElement for eccentricity connection for ETABS files - length = np.linalg.norm(np.array(geometry[1]) - np.array(geometry[0])) - for c in connections: - if c["eccentricity"]: - if ( - np.linalg.norm(np.array(c["eccentricity"]["pointOnElement"])) - > length + self.tol - ): - print( - f"{np.linalg.norm(np.array(c['eccentricity']['pointOnElement']))} > {length}" - ) - self.warnings.append( - f"Eccentricity in {item.is_a()}|{item.id()} corrected" - ) - c["eccentricity"]["pointOnElement"][0] = length - # End <-- - if transformation: - geometry = self.transform_vectors(geometry, transformation) - orientation = self.transform_vectors( - orientation, transformation, include_translation=False - ) - for c in connections: - c["orientation"] = self.transform_vectors( - c["orientation"], transformation, include_translation=False - ) - if c["eccentricity"]: - c["eccentricity"]["vector"] = self.transform_vectors( - c["eccentricity"]["vector"], - transformation, - include_translation=False, - ) - - return { - "referenceName": f"{item.is_a()}|{item.id()}", - "name": item.Name, - "id": item.GlobalId, - "geometryType": "line", - "predefinedType": item.PredefinedType, - "geometry": geometry, - "orientation": orientation, - "material": materialId, - "profile": profileId, - "connections": connections, - } - - elif item.is_a("IfcStructuralSurfaceMember"): - representation = self.get_representation(item, "Face") - material = self.get_material_profile(item) - if not representation: - self.warnings.append( - f"No representation defined for {item.is_a()}|{item.id()}. Member excluded" - ) - return - if not material: - self.warnings.append(f"No material defined for {item.is_a()}|{item.id()}") - materialId = None - else: - materialId = material.is_a() + "|" + str(material.id()) - - geometry = self.get_geometry(representation) - orientation = self.get_2D_orientation(representation) - connections = self.get_connection_data(item.ConnectedBy) - for conn in connections: - if not conn["orientation"]: - conn["orientation"] = orientation - if transformation: - geometry = self.transform_vectors(geometry, transformation) - orientation = self.transform_vectors( - orientation, transformation, include_translation=False - ) - for c in connections: - c["orientation"] = self.transform_vectors( - c["orientation"], transformation, include_translation=False - ) - - return { - "referenceName": f"{item.is_a()}|{item.id()}", - "name": item.Name, - "id": item.GlobalId, - "geometryType": "surface", - "predefinedType": item.PredefinedType, - "thickness": item.Thickness, - "geometry": geometry, - "orientation": orientation, - "material": materialId, - "connections": connections, - } - - elif item.is_a("IfcStructuralPointConnection"): - representation = self.get_representation(item, "Vertex") - if not representation: - self.warnings.append( - f"No representation defined for {item.is_a()}|{item.id()}. Member excluded" - ) - return - - geometry = self.get_geometry(representation) - orientation = self.get_0D_orientation(item.ConditionCoordinateSystem) - if not orientation: - orientation = np.eye(3).tolist() - if transformation: - geometry = self.transform_vectors(geometry, transformation) - orientation = self.transform_vectors( - orientation, transformation, include_translation=False - ) - - return { - "referenceName": f"{item.is_a()}|{item.id()}", - "name": item.Name, - "id": item.GlobalId, - "geometryType": "point", - "geometry": geometry, - "orientation": orientation, - "appliedCondition": self.get_connection_input(item, "point"), - "relatedElements": [ - f"{con.is_a()}|{con.id()}" for con in item.ConnectsStructuralMembers - ], - } - - elif item.is_a("IfcStructuralCurveConnection"): - representation = self.get_representation(item, "Edge") - if not representation: - self.warnings.append( - f"No representation defined for {item.is_a()}|{item.id()}. Member excluded" - ) - return - - geometry = self.get_geometry(representation) - orientation = self.get_1D_orientation(geometry, item.Axis) - if not orientation: - orientation = np.eye(3).tolist() - if transformation: - geometry = self.transform_vectors(geometry, transformation) - orientation = self.transform_vectors( - orientation, transformation, include_translation=False - ) - - return { - "referenceName": f"{item.is_a()}|{item.id()}", - "name": item.Name, - "id": item.GlobalId, - "geometryType": "line", - "geometry": geometry, - "orientation": orientation, - "appliedCondition": self.get_connection_input(item, "line"), - "relatedElements": [ - f"{con.is_a()}|{con.id()}" for con in item.ConnectsStructuralMembers - ], - } - - def get_transformation(self, placement): - if not placement: - return None - if placement.is_a("IfcLocalPlacement"): - if placement.PlacementRelTo: - print( - "Warning! Object Placement with PlacementRelTo attribute is not supported and will be neglected" - ) - axes = placement.RelativePlacement - location = np.array(self.get_coordinate(axes.Location)) - if axes.Axis and axes.RefDirection: - xAxis = np.array( - axes.RefDirection.DirectionRatios - ) # this can be not accurate (in the xz plane) - zAxis = np.array(axes.Axis.DirectionRatios) - zAxis /= np.linalg.norm(zAxis) - yAxis = np.cross(zAxis, xAxis) - yAxis /= np.linalg.norm(yAxis) - xAxis = np.cross(yAxis, zAxis) - xAxis /= np.linalg.norm(xAxis) - else: - if np.allclose(location, np.array([0.0, 0.0, 0.0])): - return None - xAxis = np.array([1.0, 0.0, 0.0]) - yAxis = np.array([0.0, 1.0, 0.0]) - zAxis = np.array([0.0, 0.0, 1.0]) - if ( - np.allclose(location, np.array([0.0, 0.0, 0.0])) - and np.allclose(xAxis, np.array([1.0, 0.0, 0.0])) - and np.allclose(yAxis, np.array([0.0, 1.0, 0.0])) - and np.allclose(zAxis, np.array([0.0, 0.0, 1.0])) - ): - return None - return { - "location": location, - "rotationMatrix": np.array([xAxis, yAxis, zAxis]).transpose(), - } - else: - print( - f"Warning! Object Placement is of type {placement.is_a()}, which is not supported. Default considered" - ) - return None - - def get_representation(self, element, rep_type): - if not element.Representation: - return None - for representation in element.Representation.Representations: - rep = self.get_specific_representation(representation, "Reference", rep_type) - if rep: - return rep - else: - # print("Trying without rep identifier") - for representation in element.Representation.Representations: - rep = self.get_specific_representation(representation, None, rep_type) - if rep: - return rep - - def get_specific_representation(self, representation, rep_id, rep_type): - if ( - representation.RepresentationIdentifier == rep_id or rep_id is None - ) and representation.RepresentationType == rep_type: - return representation - if representation.RepresentationType == "MappedRepresentation": - return self.get_specific_representation( - representation.Items[0].MappingSource.MappedRepresentation, - rep_id, - rep_type, - ) - - def get_geometry(self, representation): - # Maybe IfcOpenShell can use create_shape here to simplify this, but - # supposedly structural models are very simple anyway, so perhaps we - # can do without it. - item = representation.Items[0] - if item.is_a("IfcEdge"): + def get_elements(self, model: ios.entity_instance | None = None): + if model is not None: return [ - self.get_coordinate(item.EdgeStart.VertexGeometry), - self.get_coordinate(item.EdgeEnd.VertexGeometry), + item for item in ifcopenshell.util.element.get_grouped_by(model) if item.is_a("IfcStructuralMember") ] + return self.file.by_type("IfcStructuralMember") - elif item.is_a("IfcFaceSurface"): - edges = item.Bounds[0].Bound.EdgeList - coords = [] - for edge in edges: - coords.append( - self.get_coordinate(edge.EdgeElement.EdgeStart.VertexGeometry) - ) - return coords + def get_connections(self, model: ios.entity_instance | None = None): + if model is not None: + return [ + item for item in ifcopenshell.util.element.get_grouped_by(model) if item.is_a("IfcStructuralConnection") + ] + return self.file.by_type("IfcStructuralConnection") - elif item.is_a("IfcVertexPoint"): - return self.get_coordinate(item.VertexGeometry) + ## functions related to loads and load groups + def get_actions( + self, + load_group: ios.entity_instance | None = None, + element: ios.entity_instance | None = None, + load_group_ids: List[int] | None = None, + ): + if load_group is not None: + actions = [ + item + for item in ifcopenshell.util.element.get_grouped_by(load_group) + if item.is_a("IfcStructuralAction") + ] + return actions - def get_coordinate(self, point): - if point.is_a("IfcCartesianPoint"): - return list(point.Coordinates) + if element is not None: + actions = [ + item.RelatedStructuralActivity + for item in element.AssignedStructuralActivity + if item.RelatedStructuralActivity.is_a("IfcStructuralAction") + ] + if load_group_ids is None: + return actions - def get_0D_orientation(self, axes): - if axes and axes.Axis and axes.RefDirection: - xAxis = np.array( - axes.RefDirection.DirectionRatios - ) # this can be not strictly perpendicular (in the xz plane) - zAxis = np.array(axes.Axis.DirectionRatios) - zAxis /= np.linalg.norm(zAxis) - yAxis = np.cross(zAxis, xAxis) - yAxis /= np.linalg.norm(yAxis) - xAxis = np.cross(yAxis, zAxis) - xAxis /= np.linalg.norm(xAxis) + actions = [action for action in actions if self.get_load_group(action).id() in load_group_ids] + return actions - return [xAxis.tolist(), yAxis.tolist(), zAxis.tolist()] - else: # return None and copy the elements orientation + return [item for item in self.file.by_type("IfcStructuralAction")] + + def get_load_group(self, action: ios.entity_instance): + if action.HasAssignments is None: return None + load_groups = [ + item.RelatingGroup + for item in action.HasAssignments + if (item.RelatingGroup.is_a("IfcStructuralLoadGroup") and item.RelatingGroup.PredefinedType == "LOAD_GROUP") + ] + if len(load_groups): + return load_groups[0] + return None - def get_1D_orientation(self, geometry, zAxis): - xAxis = np.array(geometry[1]) - np.array(geometry[0]) - xAxis /= np.linalg.norm(xAxis) - zAxis = np.array( - zAxis.DirectionRatios - ) # this can be not strictly perpendicular (in the xz plane) - yAxis = np.cross(zAxis, xAxis) - yAxis /= np.linalg.norm(yAxis) - zAxis = np.cross(xAxis, yAxis) - zAxis /= np.linalg.norm(zAxis) + def get_reactions( + self, result_group: ios.entity_instance | None = None, element: ios.entity_instance | None = None + ): + if result_group is not None: + reactions = [ + item + for item in ifcopenshell.util.element.get_grouped_by(result_group) + if item.is_a("IfcStructuralReaction") + ] + return reactions - return [xAxis.tolist(), yAxis.tolist(), zAxis.tolist()] + if element is not None: + reactions = [ + item.RelatedStructuralActivity + for item in element.AssignedStructuralActivity + if item.RelatedStructuralActivity.is_a("IfcStructuralReaction") + ] + return reactions - def get_2D_orientation(self, representation): - item = representation.Items[0] - if item.is_a("IfcFaceSurface"): - axes = item.FaceSurface.Position - orientation = self.get_0D_orientation(axes) - if not orientation: - self.warnings.append( - f"No local placement for Plane related to {item.is_a()}|{item.id()}. A unit orientation is considered" + return [item for item in self.file.by_type("IfcStructuralReaction")] + + def get_load_groups(self, load_case: ios.entity_instance | None = None, model: ios.entity_instance | None = None): + if load_case is not None: + load_groups = [ + item + for item in ifcopenshell.util.element.get_grouped_by(load_case) + if (item.is_a("IfcStructuralLoadGroup") and item.PredefinedType == "LOAD_GROUP") + ] + return load_groups + + if model is not None: + load_group_set = set() + for load_case in self.get_load_cases(model): + load_groups = set( + [ + item + for item in ifcopenshell.util.element.get_grouped_by(load_case) + if (item.is_a("IfcStructuralLoadGroup") and item.PredefinedType == "LOAD_GROUP") + ] ) - return np.eye(3).tolist() - if not item.SameSense: - orientation = [[-v for v in vec] for vec in orientation] - return orientation + load_group_set = load_group_set.union(load_groups) + load_groups = list(load_group_set) + load_groups.sort(key=lambda x: x.id()) - def transform_vectors(self, geometry, trsf, include_translation=True): - if not any( - isinstance(el, list) for el in geometry - ): # single point which contains no list - geometry = [geometry] - globalGeometry = [] + return load_groups - for p in geometry: - gp = trsf["rotationMatrix"].dot(np.array(p)) - if include_translation: - gp += trsf["location"] - globalGeometry.append(gp.tolist()) + return [item for item in self.file.by_type("IfcStructuralLoadGroup") if item.PredefinedType == "LOAD_GROUP"] - if len(globalGeometry) == 1: # single point - globalGeometry = globalGeometry[0] + ## functions related to load cases and load combinations + def get_load_cases(self, model: ios.entity_instance | None = None, load_group: ios.entity_instance | None = None): + if model is not None: + load_case_set = set(self.get_analysis_load_cases(model)) + for comb in self.get_load_combinations(model): + load_cases = set( + [ + item + for item in ifcopenshell.util.element.get_grouped_by(comb) + if (item.is_a("IfcStructuralLoadCase") and item.PredefinedType == "LOAD_CASE") + ] + ) + load_case_set = load_case_set.union(load_cases) + load_cases = list(load_case_set) + load_cases.sort(key=lambda x: x.id()) - return globalGeometry + return load_cases - def get_material_profile(self, element): - if not element.HasAssociations: - return None - for association in element.HasAssociations: - if not association.is_a("IfcRelAssociatesMaterial"): - continue - material = association.RelatingMaterial - if material.is_a("IfcMaterialProfileSet"): - # For now, we only deal with a single profile - return material.MaterialProfiles[0] - if material.is_a("IfcMaterialProfileSetUsage"): - return material.ForProfileSet.MaterialProfiles[0] - if material.is_a("IfcMaterial"): - return material + if load_group is not None: + if load_group.HasAssignments is None: + return [] + load_cases = [ + item.RelatingGroup + for item in load_group.HasAssignments + if ( + item.RelatingGroup.is_a("IfcStructuralLoadCase") + and item.RelatingGroup.PredefinedType == "LOAD_CASE" + ) + ] + return load_cases - def get_material_properties(self, material): - psets = material.HasProperties + return [item for item in self.file.by_type("IfcStructuralLoadCase") if item.PredefinedType == "LOAD_CASE"] - if self.get_pset_properties(psets, "Pset_MaterialMechanical"): - mechProps = self.get_pset_properties(psets, "Pset_MaterialMechanical") - else: - mechProps = self.get_pset_properties(psets, None) + def get_load_combinations(self, model: ios.entity_instance | None = None): + if model is not None: + if model.LoadedBy is None: + return [] + load_groups = [ + item + for item in model.LoadedBy + if (item.is_a("IfcStructuralLoadGroup") and item.PredefinedType == "LOAD_COMBINATION") + ] + load_groups.sort(key=lambda x: x.id()) + return load_groups - if self.get_pset_properties(psets, "Pset_MaterialCommon"): - commonProps = self.get_pset_properties(psets, "Pset_MaterialCommon") - else: - commonProps = self.get_pset_properties(psets, None) - - return { - "referenceName": material.is_a() + "|" + str(material.id()), - "name": material.Name, - "category": material.Category, - "mechProps": mechProps, - "commonProps": commonProps, - } - - def get_pset_property(self, psets, pset_name, prop_name): - for pset in psets: - if pset.Name == pset_name or pset_name is None: - for prop in pset.Properties: - if prop.Name == prop_name: - return prop.NominalValue.wrappedValue - - def get_pset_properties(self, psets, pset_name): - for pset in psets: - if pset.Name == pset_name or pset_name is None: - d = {} - for prop in pset.Properties: - propName = prop.Name[0].lower() + prop.Name[1:] - d[propName] = prop.NominalValue.wrappedValue - return d - - def get_profile_properties(self, profile): - if profile.is_a("IfcRectangleProfileDef"): - return { - "referenceName": profile.is_a() + "|" + str(profile.id()), - "profileName": profile.ProfileName, - "profileType": profile.ProfileType, - "profileShape": "rectangular", - "xDim": profile.XDim, - "yDim": profile.YDim, - } - - if profile.is_a("IfcIShapeProfileDef"): - psets = profile.HasProperties - - if self.get_pset_properties(psets, "Pset_ProfileMechanical"): - mechProps = self.get_pset_properties(psets, "Pset_ProfileMechanical") - else: - mechProps = self.get_i_section_properties(profile, "iSymmetrical") - - return { - "referenceName": f"{profile.is_a()}|{profile.id()}", - "profileName": profile.ProfileName, - "profileType": profile.ProfileType, - "profileShape": "iSymmetrical", - "mechProps": mechProps, - "commonProps": { - "flangeThickness": profile.FlangeThickness, - "webThickness": profile.WebThickness, - "overallDepth": profile.OverallDepth, - "overallWidth": profile.OverallWidth, - "filletRadius": profile.FilletRadius, - }, - } - - def get_connection_data(self, itemList): return [ - { - "referenceName": f"{rel.is_a()}|{rel.id()}", - "id": rel.GlobalId, - "relatingElement": f"{rel.RelatingStructuralMember.is_a()}|{rel.RelatingStructuralMember.id()}", - "relatedConnection": f"{rel.RelatedStructuralConnection.is_a()}|{rel.RelatedStructuralConnection.id()}", - "orientation": self.get_0D_orientation(rel.ConditionCoordinateSystem), - "appliedCondition": self.get_connection_input( - rel, - self.get_geometry_type_from_connection( - rel.RelatedStructuralConnection - ), - ), - "eccentricity": None - if not rel.is_a("IfcRelConnectsWithEccentricity") - else { - "vector": [ - 0.0 - if not rel.ConnectionConstraint.EccentricityInX - else rel.ConnectionConstraint.EccentricityInX, - 0.0 - if not rel.ConnectionConstraint.EccentricityInY - else rel.ConnectionConstraint.EccentricityInY, - 0.0 - if not rel.ConnectionConstraint.EccentricityInZ - else rel.ConnectionConstraint.EccentricityInZ, - ], - "pointOnElement": self.get_coordinate( - rel.ConnectionConstraint.PointOnRelatingElement - ), - }, - } - for rel in itemList + item for item in self.file.by_type("IfcStructuralLoadGroup") if item.PredefinedType == "LOAD_COMBINATION" ] - def get_geometry_type_from_connection(self, connection): - if connection.is_a("IfcStructuralPointConnection"): - return "point" - if connection.is_a("IfcStructuralCurveConnection"): - return "line" - if connection.is_a("IfcStructuralSurfaceConnection"): - return "surface" + def get_combination_assignments(self, load_combination: ios.entity_instance): + return [ + rel for rel in self.file.by_type("IfcRelAssignsToGroup") if rel.RelatingGroup.id() == load_combination.id() + ] - def get_connection_input(self, connection, geometryType): - if connection.AppliedCondition: - if geometryType == "point": - return { - "dx": connection.AppliedCondition.TranslationalStiffnessX.wrappedValue, - "dy": connection.AppliedCondition.TranslationalStiffnessY.wrappedValue, - "dz": connection.AppliedCondition.TranslationalStiffnessZ.wrappedValue, - "drx": connection.AppliedCondition.RotationalStiffnessX.wrappedValue, - "dry": connection.AppliedCondition.RotationalStiffnessY.wrappedValue, - "drz": connection.AppliedCondition.RotationalStiffnessZ.wrappedValue, + def get_analysis_load_cases(self, model: ios.entity_instance): + if model.LoadedBy is None: + return [] + return [ + item + for item in model.LoadedBy + if (item.is_a("IfcStructuralLoadCase") and item.PredefinedType == "LOAD_CASE") + ] + + def get_analysis_cases(self, model: ios.entity_instance): + if model.LoadedBy is None: + return [] + return list(model.LoadedBy) + + def get_result_cases(self, model: ios.entity_instance): + if model.HasResults is None: + return [] + return list(model.HasResults) + + # expose GET/PARSE api functions + # functions to convert to json + def get_ref_id(self, entity: ios.entity_instance): + return f"{entity.is_a()}|{entity.id()}" + + def parse_transformation_matrix(self, matrix: np.ndarray): + x_axis = np.round([v[0] for v in matrix][:3], 6) + y_axis = np.round([v[1] for v in matrix][:3], 6) + z_axis = np.round([v[2] for v in matrix][:3], 6) + + origin = np.round([v[3] for v in matrix][:3], 6) + orientation = np.array([x_axis, y_axis, z_axis]) + + return origin, orientation + + def parse_representation(self, representation: ios.entity_instance): + repr_item = representation.Items[0] + if representation.RepresentationType == "Vertex": + geometry = repr_item.VertexGeometry.Coordinates + + elif representation.RepresentationType == "Edge": + geometry = [ + repr_item.EdgeStart.VertexGeometry.Coordinates, + repr_item.EdgeEnd.VertexGeometry.Coordinates, + ] + + elif representation.RepresentationType == "Face": + geometry = [x.EdgeStart.VertexGeometry.Coordinates for x in repr_item.Bounds[0].Bound.EdgeList] + + else: + print(representation) + return geometry + + def parse_material(self, material: ios.entity_instance): + data = {k: v for k, v in material.get_info().items() if k in self.material_keys} + data["ref_id"] = self.get_ref_id(material) + + psets = ifcopenshell.util.element.get_psets(material) + properties = dict() + for _, pset_properties in psets.items(): + properties |= pset_properties + if len(properties): + data["properties"] = {k: v for k, v in properties.items() if k in self.material_property_keys} + if "PoissonRatio" not in data["properties"]: + if "ShearModulus" in data["properties"]: + PoissonRatio = round( + data["properties"]["YoungModulus"] / 2.0 / data["properties"]["ShearModulus"] - 1, 3 + ) + else: + PoissonRatio = 0.0 + data["properties"]["PoissonRatio"] = PoissonRatio + else: + data["properties"] = None + + return data + + def parse_profile(self, profile: ios.entity_instance): + data = profile.get_info() + data["ref_id"] = self.get_ref_id(profile) + + psets = ifcopenshell.util.element.get_psets(profile) + properties = dict() + for _, pset_properties in psets.items(): + properties |= pset_properties + if len(properties): + data["properties"] = {k: v for k, v in properties.items() if k in self.profile_property_keys} + else: + if profile.is_a("IfcIShapeProfileDef"): + tf = profile.FlangeThickness + tw = profile.WebThickness + h = profile.OverallDepth + b = profile.OverallWidth + + A = b * h - (b - tw) * (h - 2 * tf) + Iy = b * (h**3) / 12 - (b - tw) * ((h - 2 * tf) ** 3) / 12 + Iz = (2 * tf) * (b**3) / 12 + (h - 2 * tf) * (tw**3) / 12 + Jx = 1 / 3 * ((h - tf) * (tw**3) + 2 * b * (tf**3)) + + data["properties"] = { + "CrossSectionArea": A, + "MomentOfInertiaY": Iy, + "MomentOfInertiaZ": Iz, + "TorsionalConstantX": Jx, } - if geometryType == "line": - return { - "dx": connection.AppliedCondition.TranslationalStiffnessByLengthX.wrappedValue, - "dy": connection.AppliedCondition.TranslationalStiffnessByLengthY.wrappedValue, - "dz": connection.AppliedCondition.TranslationalStiffnessByLengthZ.wrappedValue, - "drx": connection.AppliedCondition.RotationalStiffnessByLengthX.wrappedValue, - "dry": connection.AppliedCondition.RotationalStiffnessByLengthY.wrappedValue, - "drz": connection.AppliedCondition.RotationalStiffnessByLengthZ.wrappedValue, - } + else: + data["properties"] = None - if geometryType == "surface": - return { - "dx": connection.AppliedCondition.TranslationalStiffnessByAreaX.wrappedValue, - "dy": connection.AppliedCondition.TranslationalStiffnessByAreaY.wrappedValue, - "dz": connection.AppliedCondition.TranslationalStiffnessByAreaZ.wrappedValue, - } + return data - return connection.AppliedCondition - - def get_i_section_properties(self, profile, profileShape): - if profileShape == "iSymmetrical": - tf = profile.FlangeThickness - tw = profile.WebThickness - h = profile.OverallDepth - b = profile.OverallWidth - - A = b * h - (b - tw) * (h - 2 * tf) - Iy = b * (h ** 3) / 12 - (b - tw) * ((h - 2 * tf) ** 3) / 12 - Iz = (2 * tf) * (b ** 3) / 12 + (h - 2 * tf) * (tw ** 3) / 12 - Jx = 1 / 3 * ((h - tf) * (tw ** 3) + 2 * b * (tf ** 3)) + def parse_applied_condition( + self, condition: ios.entity_instance, geometry_type: str, should_fix_condition: bool = False + ): + if condition is None: + if should_fix_condition: + return self.fixed_conditions[geometry_type] + else: + return None + if geometry_type == "Vertex": return { - "crossSectionArea": A, - "momentOfInertiaY": Iy, - "momentOfInertiaZ": Iz, - "torsionalConstantX": Jx, + "dx": condition.TranslationalStiffnessX.wrappedValue, + "dy": condition.TranslationalStiffnessY.wrappedValue, + "dz": condition.TranslationalStiffnessZ.wrappedValue, + "drx": condition.RotationalStiffnessX.wrappedValue, + "dry": condition.RotationalStiffnessY.wrappedValue, + "drz": condition.RotationalStiffnessZ.wrappedValue, } + if geometry_type == "Edge": + return { + "dx": condition.TranslationalStiffnessByLengthX.wrappedValue, + "dy": condition.TranslationalStiffnessByLengthY.wrappedValue, + "dz": condition.TranslationalStiffnessByLengthZ.wrappedValue, + "drx": condition.RotationalStiffnessByLengthX.wrappedValue, + "dry": condition.RotationalStiffnessByLengthY.wrappedValue, + "drz": condition.RotationalStiffnessByLengthZ.wrappedValue, + } -if __name__ == "__main__": - fileNames = [ - "cantilever_01", - "portal_01", - "grid_of_beams", - "slab_01", - "structure_01", - "building_02", - ] - files = fileNames + if geometry_type == "Face": + return { + "dx": condition.TranslationalStiffnessByAreaX.wrappedValue, + "dy": condition.TranslationalStiffnessByAreaY.wrappedValue, + "dz": condition.TranslationalStiffnessByAreaZ.wrappedValue, + } - for fileName in files: - BASE_PATH = Path( - "/home/jesusbill/Dev-Projects/github.com/IfcOpenShell/analysis-models/ifcFiles/" + def parse_element(self, element: ios.entity_instance): + data = {k: v for k, v in element.get_info().items() if k in self.member_keys} + + data["ObjectType"] = ifcopenshell.util.element.get_predefined_type(element) + data["ref_id"] = self.get_ref_id(element) + representation = ifcopenshell.util.representation.get_representation(element, self.get_context()) + repr_item = representation.Items[0] + data["geometry_type"] = representation.RepresentationType + data["geometry"] = self.parse_representation(representation) + + if element.is_a("IfcStructuralCurveMember"): + placement = ifcopenshell.util.placement.a2p( + data["geometry"][0], + element.Axis.DirectionRatios, + [c2 - c1 for c1, c2 in zip(data["geometry"][0], data["geometry"][1])], + ) + + elif element.is_a("IfcStructuralSurfaceMember"): + placement = ifcopenshell.util.placement.get_axis2placement(repr_item.FaceSurface.Position) + + origin, orientation = self.parse_transformation_matrix(placement) + data["origin"] = origin + data["orientation"] = orientation + if element.is_a("IfcStructuralSurfaceMember") and not repr_item.SameSense: + data["orientation"] *= -1 + + materialset = ifcopenshell.util.element.get_material(element, should_skip_usage=True) + if element.is_a() == "IfcStructuralCurveMember": + data["material"] = materialset.MaterialProfiles[0].Material + data["profile"] = materialset.MaterialProfiles[0].Profile + # data["profile"] = ifcopenshell.util.shape.get_profiles(element)[0] + elif element.is_a() == "IfcStructuralSurfaceMember": + if materialset.is_a() == "IfcMaterial": + data["material"] = materialset + else: + data["material"] = materialset.MaterialLayers[0].Material + + return data + + def parse_connection(self, connection: ios.entity_instance): + data = {k: v for k, v in connection.get_info().items() if k in self.connection_keys} + + data["ObjectType"] = ifcopenshell.util.element.get_predefined_type(connection) + data["ref_id"] = self.get_ref_id(connection) + representation = ifcopenshell.util.representation.get_representation(connection, self.get_context()) + repr_item = representation.Items[0] + data["geometry_type"] = representation.RepresentationType + data["geometry"] = self.parse_representation(representation) + + if connection.is_a("IfcStructuralPointConnection"): + if connection.ConditionCoordinateSystem is not None: + placement = ifcopenshell.util.placement.get_axis2placement(connection.ConditionCoordinateSystem) + else: + placement = np.eye(4) + for i, v in enumerate(placement[:3]): + v[3] = data["geometry"][i] + + if connection.is_a("IfcStructuralCurveConnection"): + placement = ifcopenshell.util.placement.a2p( + data["geometry"][0], + connection.Axis.DirectionRatios, + [c2 - c1 for c1, c2 in zip(data["geometry"][0], data["geometry"][1])], + ) + + elif connection.is_a("IfcStructuralSurfaceConnection"): + placement = ifcopenshell.util.placement.get_axis2placement(repr_item.FaceSurface.Position) + + origin, orientation = self.parse_transformation_matrix(placement) + data["origin"] = origin + data["orientation"] = orientation + if connection.is_a("IfcStructuralSurfaceConnection") and not repr_item.SameSense: + data["orientation"] *= -1 + + data["appliedCondition"] = self.parse_applied_condition(connection.AppliedCondition, data["geometry_type"]) + return data + + def parse_element_connections(self, element: Dict, connections: List[Dict]): + ifc_element = self.file.by_id(element["id"]) + connection_ids = [c["id"] for c in connections] + rels = [rel for rel in ifc_element.ConnectedBy if rel.RelatedStructuralConnection.id() in connection_ids] + + data = [dict() for _ in rels] + for i, rel in enumerate(rels): + data[i]["ref_id"] = self.get_ref_id(rel) + data[i]["related_connection"] = self.get_ref_id(rel.RelatedStructuralConnection) + data[i]["relating_element"] = self.get_ref_id(rel.RelatingStructuralMember) + + if rel.ConditionCoordinateSystem is None: + data[i]["orientation"] = element["orientation"] + else: + placement = ifcopenshell.util.placement.get_axis2placement(rel.ConditionCoordinateSystem) + _, orientation = self.parse_transformation_matrix(placement) + data[i]["orientation"] = element["orientation"].dot(orientation) + + conn_repr = ifcopenshell.util.representation.get_representation( + rel.RelatedStructuralConnection, self.get_context() + ) + data[i]["geometry_type"] = conn_repr.RepresentationType + data[i]["appliedCondition"] = self.parse_applied_condition( + rel.AppliedCondition, data[i]["geometry_type"], should_fix_condition=True + ) + + if rel.is_a("IfcRelConnectsWithEccentricity") and rel.RelatedStructuralConnection.is_a( + "IfcStructuralPointConnection" + ): + point_on_element = rel.ConnectionConstraint.PointOnRelatingElement.Coordinates + element_length = round( + np.linalg.norm(np.array(element["geometry"][0]) - np.array(element["geometry"][1])), 6 + ) + x_local = min(max(0.0, point_on_element[0]), element_length) + point_on_element = ( + np.array(element["geometry"][0]) + + (np.array(element["geometry"][1]) - np.array(element["geometry"][0])) * x_local / element_length + ) + point_on_element = np.round(point_on_element, 6) + data[i]["eccentricity"] = {"point_on_element": point_on_element} + else: + data[i]["eccentricity"] = ( + rel.ConnectionConstraint.get_info(recursive=True) + if rel.is_a("IfcRelConnectsWithEccentricity") + else None + ) + + return data + + def parse_element_loads( + self, + element: Dict, + load_group_ids: List[int], + load_cases: List[ios.entity_instance], + ): + ifc_element = self.file.by_id(element["id"]) + actions = self.get_actions(element=ifc_element, load_group_ids=load_group_ids) + if not len(actions): + return None + + if element["geometry_type"] in ["Vertex", "Edge"]: + data = { + "actions": [], + "loadGroups": [], + "loadsLC": { + "FX": np.array([0.0 for _ in load_cases]), + "FY": np.array([0.0 for _ in load_cases]), + "FZ": np.array([0.0 for _ in load_cases]), + "MX": np.array([0.0 for _ in load_cases]), + "MY": np.array([0.0 for _ in load_cases]), + "MZ": np.array([0.0 for _ in load_cases]), + }, + "loadsCOMB": { + "FX": None, + "FY": None, + "FZ": None, + "MX": None, + "MY": None, + "MZ": None, + }, + } + elif element["geometry_type"] == "Face": + data = { + "actions": [], + "loadGroups": [], + "loadsLC": { + "FX": np.array([0.0 for _ in load_cases]), + "FY": np.array([0.0 for _ in load_cases]), + "FZ": np.array([0.0 for _ in load_cases]), + }, + "loadsCOMB": { + "FX": None, + "FY": None, + "FZ": None, + }, + } + + for action in actions: + self.add_action_loads(element, action, data, load_cases) + + loadsLC = deepcopy(data["loadsLC"]) + loadsCOMB = deepcopy(data["loadsCOMB"]) + for key, load in data["loadsLC"].items(): + if np.max(np.abs(load)) == 0.0: + del loadsLC[key] + del loadsCOMB[key] + if len(loadsLC) == 0: + print( + f"Actions {[self.get_ref_id(action) for action in actions]} for element {element['ref_id']} not accounted for in the load assignments" + ) + return None + else: + data["loadsLC"] = loadsLC + data["loadsCOMB"] = loadsCOMB + + return data + + def add_action_loads( + self, + element: Dict, + action: ios.entity_instance, + data: Dict, + load_cases: List[ios.entity_instance], + ): + load_group = self.get_load_group(action) + load_group_coeff = 1.0 if load_group.Coefficient is None else load_group.Coefficient + active_load_case_ids = [lc.id() for lc in self.get_load_cases(load_group=load_group)] + load = action.AppliedLoad + + data["actions"].append(action.get_info() | {"AppliedLoad": action.AppliedLoad.get_info()}) + if element["geometry_type"] in ["Vertex", "Edge"]: + if action.is_a("IfcStructuralPointAction") and load.is_a("IfcStructuralLoadSingleForce"): + FX = tempFX = load.ForceX if load.ForceX is not None else 0.0 + FY = tempFY = load.ForceY if load.ForceY is not None else 0.0 + FZ = tempFZ = load.ForceZ if load.ForceZ is not None else 0.0 + MX = tempMX = load.MomentX if load.MomentX is not None else 0.0 + MY = tempMY = load.MomentY if load.MomentY is not None else 0.0 + MZ = tempMZ = load.MomentZ if load.MomentZ is not None else 0.0 + if action.GlobalOrLocal == "LOCAL_COORDS": + # adjust for global/local + FX = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 0]) + FY = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 1]) + FZ = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 2]) + MX = np.array([tempMX, tempMY, tempMZ]).dot(element["orientation"][:, 0]) + MY = np.array([tempMX, tempMY, tempMZ]).dot(element["orientation"][:, 1]) + MZ = np.array([tempMX, tempMY, tempMZ]).dot(element["orientation"][:, 2]) + + elif action.is_a("IfcStructuralLinearAction") and load.is_a("IfcStructuralLoadLinearForce"): + FX = tempFX = load.LinearForceX if load.LinearForceX is not None else 0.0 + FY = tempFY = load.LinearForceY if load.LinearForceY is not None else 0.0 + FZ = tempFZ = load.LinearForceZ if load.LinearForceZ is not None else 0.0 + MX = tempMX = load.LinearMomentX if load.LinearMomentX is not None else 0.0 + MY = tempMY = load.LinearMomentY if load.LinearMomentY is not None else 0.0 + MZ = tempMZ = load.LinearMomentZ if load.LinearMomentZ is not None else 0.0 + if action.GlobalOrLocal == "LOCAL_COORDS": + # adjust for global/local and projected/true + FX = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 0]) + FY = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 1]) + FZ = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 2]) + MX = np.array([tempMX, tempMY, tempMZ]).dot(element["orientation"][:, 0]) + MY = np.array([tempMX, tempMY, tempMZ]).dot(element["orientation"][:, 1]) + MZ = np.array([tempMX, tempMY, tempMZ]).dot(element["orientation"][:, 2]) + force_projection_coeff = 1.0 + moment_projection_coeff = 1.0 + else: + if action.ProjectedOrTrue == "PROJECTED_LENGTH": + if force_length := np.linalg.norm(np.array([FX, FY, FZ])): + unit_force = np.array([FX, FY, FZ]) / force_length + + force_projection_coeff = 1 - round(abs(unit_force.dot(element["orientation"][0])), 4) + assert force_projection_coeff >= 0, f"{force_projection_coeff}" + else: + force_projection_coeff = 1.0 + + if moment_length := np.linalg.norm(np.array([MX, MY, MZ])): + unit_moment = np.array([MX, MY, MZ]) / moment_length + moment_projection_coeff = 1 - round(abs(unit_moment.dot(element["orientation"][0])), 4) + assert moment_projection_coeff >= 0, f"{moment_projection_coeff}" + else: + moment_projection_coeff = 1.0 + else: + force_projection_coeff = 1.0 + moment_projection_coeff = 1.0 + + for iLC, load_case in enumerate(load_cases): + if load_case.id() in active_load_case_ids: + load_case_coeff = 1.0 if load_case.Coefficient is None else load_case.Coefficient + data["loadGroups"].append(load_group.Name) + data["loadsLC"]["FX"][iLC] += FX * load_group_coeff * load_case_coeff * force_projection_coeff + data["loadsLC"]["FY"][iLC] += FY * load_group_coeff * load_case_coeff * force_projection_coeff + data["loadsLC"]["FZ"][iLC] += FZ * load_group_coeff * load_case_coeff * force_projection_coeff + data["loadsLC"]["MX"][iLC] += MX * load_group_coeff * load_case_coeff * moment_projection_coeff + data["loadsLC"]["MY"][iLC] += MY * load_group_coeff * load_case_coeff * moment_projection_coeff + data["loadsLC"]["MZ"][iLC] += MZ * load_group_coeff * load_case_coeff * moment_projection_coeff + + elif element["geometry_type"] == "Face": + if action.is_a("IfcStructuralSurfaceAction") and load.is_a("IfcStructuralLoadPlanarForce"): + FX = tempFX = load.PlanarForceX if load.PlanarForceX is not None else 0.0 + FY = tempFY = load.PlanarForceY if load.PlanarForceY is not None else 0.0 + FZ = tempFZ = load.PlanarForceZ if load.PlanarForceZ is not None else 0.0 + if action.GlobalOrLocal == "LOCAL_COORDS": + # adjust for global/local and projected/true + FX = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 0]) + FY = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 1]) + FZ = np.array([tempFX, tempFY, tempFZ]).dot(element["orientation"][:, 2]) + force_projection_coeff = 1.0 + else: + if action.ProjectedOrTrue == "PROJECTED_LENGTH": + if force_length := np.linalg.norm(np.array([FX, FY, FZ])): + unit_force = np.array([FX, FY, FZ]) / force_length + force_projection_coeff = round(abs(unit_force.dot(element["orientation"][2])), 4) + assert force_projection_coeff >= 0, f"{force_projection_coeff}" + else: + force_projection_coeff = 1.0 + else: + force_projection_coeff = 1.0 + + for iLC, load_case in enumerate(load_cases): + if load_case.id() in active_load_case_ids: + load_case_coeff = 1.0 if load_case.Coefficient is None else load_case.Coefficient + data["loadGroups"].append(load_group.Name) + data["loadsLC"]["FX"][iLC] += FX * load_group_coeff * load_case_coeff * force_projection_coeff + data["loadsLC"]["FY"][iLC] += FY * load_group_coeff * load_case_coeff * force_projection_coeff + data["loadsLC"]["FZ"][iLC] += FZ * load_group_coeff * load_case_coeff * force_projection_coeff + + def parse_combination_loads(self, elements, load_combinations, load_case_ids): + for element in elements: + if element["loads"] is not None: + for key, _ in element["loads"]["loadsCOMB"].items(): + element["loads"]["loadsCOMB"][key] = np.array([0.0 for _ in load_combinations]) + + for iComb, comb in enumerate(load_combinations): + comb_factors = self.get_combination_factors(comb, load_case_ids) + # print(comb_factors) + for element in elements: + if element["loads"] is None: + continue + for key, load in element["loads"]["loadsLC"].items(): + element["loads"]["loadsCOMB"][key][iComb] = round(comb_factors.dot(load), 4) + + def get_combination_factors(self, load_combination, load_case_ids: List[int]): + comb_coeff = 1.0 if load_combination.Coefficient is None else load_combination.Coefficient + comb_factors = np.array([0.0 for _ in load_case_ids]) + for assignment in self.get_combination_assignments(load_combination): + for group in assignment.RelatedObjects: + if group.id() in load_case_ids: + iLC = load_case_ids.index(group.id()) + if assignment.is_a("IfcRelAssignsToGroupByFactor"): + comb_factors[iLC] += assignment.Factor * comb_coeff + else: + comb_factors[iLC] += 1.0 * comb_coeff + + return np.round(comb_factors, 4) + + def parse_model(self, model: ios.entity_instance): + data = {"model": {k: v for k, v in model.get_info().items() if k in self.model_keys}} + + data["model"]["ObjectType"] = ifcopenshell.util.element.get_predefined_type(model) + if model.LoadedBy: + data["model"]["LoadedBy"] = [item.Name for item in model.LoadedBy] + if model.HasResults: + data["model"]["HasResults"] = [item.Name for item in model.HasResults] + data["model"]["ref_id"] = self.get_ref_id(model) + placement = ifcopenshell.util.placement.get_local_placement(model.SharedPlacement) + origin, orientation = self.parse_transformation_matrix(placement) + data["model"]["origin"] = origin + data["model"]["orientation"] = orientation + + # elements and connections + elements = [self.parse_element(element) for element in self.get_elements(model=model)] + connections = [self.parse_connection(connection) for connection in self.get_connections(model=model)] + for element in elements: + element["connections"] = self.parse_element_connections(element, connections=connections) + + # loads, laod cases and laod combinations + load_cases = self.get_load_cases(model=model) + load_case_ids = [load_case.id() for load_case in load_cases] + load_groups = self.get_load_groups(model=model) + load_group_ids = [load_group.id() for load_group in load_groups] + for element in elements + connections: + element["loads"] = self.parse_element_loads( + element=element, + load_group_ids=load_group_ids, + load_cases=load_cases, + ) + load_combinations = self.get_load_combinations(model=model) + if len(load_combinations): + self.parse_combination_loads(elements + connections, load_combinations, load_case_ids) + + data["load_cases"] = [load_case.get_info() for load_case in load_cases] + data["load_combinations"] = [load_combination.get_info() for load_combination in load_combinations] + for combination in data["load_combinations"]: + combination["factors"] = self.get_combination_factors(self.file.by_id(combination["id"]), load_case_ids) + combination["load_cases"] = [ + item.Name + for item in ifcopenshell.util.element.get_grouped_by(self.file.by_id(combination["id"])) + if item.is_a("IfcStructuralLoadCase") + ] + + # materials and profiles + materials = set([item["material"] for item in elements]) + materialdb = dict( + zip([self.get_ref_id(mat) for mat in materials], [self.parse_material(mat) for mat in materials]) ) - ifc2ca = IFC2CA(BASE_PATH / f"{fileName}.ifc") - ifc2ca.convert() - with open(BASE_PATH / f"{fileName}.json", "w") as f: - f.write(json.dumps(ifc2ca.result, indent=4)) + for _, material in materialdb.items(): + material["related_elements"] = [ + item["ref_id"] for item in elements if material["id"] == item["material"].id() + ] + + profiles = set([item["profile"] for item in elements if item["geometry_type"] == "Edge"]) + profiledb = dict( + zip([self.get_ref_id(prof) for prof in profiles], [self.parse_profile(prof) for prof in profiles]) + ) + for _, profile in profiledb.items(): + profile["related_elements"] = [ + item["ref_id"] + for item in elements + if (item["geometry_type"] == "Edge" and profile["id"] == item["profile"].id()) + ] + data["elements"] = elements + data["connections"] = connections + data["db"] = {"materials": materialdb, "profiles": profiledb} + + return json.loads(json.dumps(data, cls=IFC2JSONEncoder)) + + # expose CRUD api functions + # functions for meshes + def get_meshes(self, model: ios.entity_instance | None = None): + if model is not None: + models = [model] + else: + models = self.get_models() + if not self.folder_path.exists(): + return [] + + meshes = [] + for model in models: + for med_file in self.folder_path.glob(f"Model_{model.id()}*.med"): + meshes.append({"name": med_file.stem, "model_id": model.id(), "med_path": str(med_file.resolve())}) + + return meshes + + def create_mesh(self, model, parameters): + model_id = model.id() + if not self.folder_path.exists(): + self.folder_path.mkdir() + + template = self.env.get_template("salome/scriptSalome.py") + mesh_name = f'Model_{model_id}_v{parameters["mesh_size"]}' + + med_path = self.folder_path / f"{mesh_name}.med" + json_path = self.folder_path / f"Model_{model_id}.json" + with json_path.open("w") as f: + json.dump(self.parse_model(model), f, indent=4) + + rendered_script = template.render( + mesh_size=parameters["mesh_size"], + json_path=str(json_path.resolve()), + med_path=str(med_path.resolve()), + mesh_name=mesh_name, + ) + + # Write the rendered script to the new location + script_path = self.folder_path / f'ScriptSalome_Model_{model_id}_v{parameters["mesh_size"]}.py' + with open(script_path, "w") as f: + f.write(rendered_script) + + if self.salome_path is not None: + executable = Path(self.salome_path) + subprocess.run(["python", executable, "-t", script_path]) + + return { + "name": mesh_name, + "model_id": model.id(), + "med_path": med_path, + } + else: + print("Salome path not provided. Mesh operation aborted.") + print(f"The salome script file path can be found in the return statement.") + return { + "name": mesh_name, + "model_id": model.id(), + "script_path": script_path, + } + + def get_run_case_labels(self, data, target): + cases = [] + if target in ["Any", "LoadCase"]: + if len(data["load_cases"]): + cases.append("LC") + + if target in ["Any", "LoadCombination"]: + if len(data["load_combinations"]): + cases.append("COMB") + + return cases + + def run_code_aster(self, mesh, target: str): + model = self.file.by_id(mesh["model_id"]) + json_path = self.folder_path / f"Model_{model.id()}.json" + + # Read data from data file + with open(json_path, "r") as f: + data = json.load(f) + cases = self.get_run_case_labels(data, target) + run_label = "_".join(cases) + + comm_path = self.folder_path / f'{mesh["name"]}_{run_label}.comm' + + constructor = CommandFileConstructor(data) + constructor.create_comm(comm_path, cases=cases) + + export_template = self.env.get_template("codeaster/export") + export_content = export_template.render( + model_name=mesh["name"], + allocated_memory=7000.0, # in MBs + time_limit=900000.0, # in seconds + cases=cases, + run_label=run_label, + ) + export_path = self.folder_path / f"export" + with open(export_path, "w") as f: + f.write(export_content) + + subprocess.run( + [ + "docker", + "run", + "-ti", + "--rm", + "-v", + f"{self.folder_path.resolve()}:/home/aster/shared", + "-w", + "/home/aster/shared", + "aethereng/ifc2ca", + "as_run", + "export", + ] + ) + return {"comm_path": comm_path, "export_path": export_path} + + # except Exception as e: + # # print(f"Exception `{e}` raised. Run operation aborted.") + # # print(f"The command and export file paths can be found in the return statement.") + # return {"Exception": e, "comm_path": comm_path, "export_path": export_path} + + def write_results_to_ifc(self, mesh, target, fields): + model = self.file.by_id(mesh["model_id"]) + json_path = self.folder_path / f"Model_{model.id()}.json" + + # Read data from data file + with open(json_path, "r") as f: + data = json.load(f) + cases = self.get_run_case_labels(data, target) + model_name = mesh["name"] + + for case_instant in cases: + rmed_path = self.folder_path / f"{model_name}_{case_instant}.rmed" + + ca2ifc.results_to_ifc(self.file, model, rmed_path, case_instant, fields, data) + + def save(self): + self.file.write(self.path) + + def save_as(self, path: os.PathLike | str): + self.file.write(path) + + def get_info(self): + models = self.get_models() + data = dict() + for model in models: + info = { + "number_of_items": len(self.get_items(model)), + "number_of_elements": len(self.get_elements(model)), + "number_of_connections": len(self.get_connections(model)), + "loads": None, + "results": None, + } + if model.LoadedBy is not None: + info["loads"] = { + "number_of_analysiscases": len(self.get_analysis_cases(model)), + "number_of_loadcases": len(self.get_analysis_load_cases(model)), + "number_of_loadcombinations": len(self.get_load_combinations(model)), + } + if model.HasResults is not None: + info["results"] = { + "number_of_resultcases": len(self.get_result_cases(model)), + } + data[self.get_ref_id(model)] = info + + return data + + def toJSON(data): + return json.dumps(data, indent=4, sort_keys=False, cls=IFC2JSONEncoder) + + +class IFC2JSONEncoder(json.JSONEncoder): + def default(self, obj): + # print(type(obj), obj) + if isinstance(obj, ios.entity_instance): + return f"{obj.is_a()}|{obj.id()}" + # return obj.get_info() + if isinstance(obj, np.ndarray): + return obj.tolist() + + # Let the base class default method raise the TypeError + return json.JSONEncoder.default(self, obj) diff --git a/src/ifc2ca/scriptCodeAster.py b/src/ifc2ca/scriptCodeAster.py index ff0ab2c272..50c8fd5ccd 100644 --- a/src/ifc2ca/scriptCodeAster.py +++ b/src/ifc2ca/scriptCodeAster.py @@ -1,5 +1,5 @@ # Ifc2CA - IFC Code_Aster utility -# Copyright (C) 2020, 2021 Ioannis P. Christovasilis +# Copyright (C) 2020, 2021, 2023, 2024 Ioannis P. Christovasilis # # This file is part of Ifc2CA. # @@ -16,20 +16,23 @@ # You should have received a copy of the GNU Lesser General Public License # along with Ifc2CA. If not, see . -import json -import numpy as np import itertools +import json from pathlib import Path +from typing import Dict + +import numpy as np +from jinja2 import Environment, FileSystemLoader flatten = itertools.chain.from_iterable -includeZeroLength1DSprings = True +includeZeroLength1DSprings = False -class COMMANDFILE: - def __init__(self, dataFilename, asterFilename): - self.dataFilename = dataFilename - self.asterFilename = asterFilename - self.create() + +class CommandFileConstructor: + def __init__(self, data: Dict): + self.data = data + self.env = Environment(loader=FileSystemLoader(Path(__file__).parent / "templates")) def getGroupName(self, name): if "|" in name: @@ -39,862 +42,6 @@ class COMMANDFILE: else: return name - def create(self): - - AccelOfGravity = 9.806 # m/sec^2 - - # Read data from input file - with open(self.dataFilename) as dataFile: - data = json.load(dataFile) - - elements = data["elements"] - connections = data["connections"] - # --> Delete this reference data and repopulate it with the objects - # while going through elements - for conn in connections: - conn["relatedElements"] = [] - self.calculateRestraints(conn) - for el in elements: - for rel in el["connections"]: - conn = [ - c for c in connections if c["referenceName"] == rel["relatedConnection"] - ][0] - rel["conn_string"] = None - if conn["geometryType"] == "point": - rel["conn_string"] = "_0DC_" - rel["springGroupName"] = ( - self.getGroupName(rel["relatingElement"]) - + "_1DS_" - + self.getGroupName(rel["relatedConnection"]) - ) - if conn["geometryType"] == "line": - rel["conn_string"] = "_1DC_" - rel["springGroupName"] = None - if conn["geometryType"] == "surface": - rel["conn_string"] = "_2DC_" - rel["springGroupName"] = None - - rel["groupName1"] = ( - self.getGroupName(rel["relatingElement"]) - + rel["conn_string"] - + self.getGroupName(rel["relatedConnection"]) - ) - if rel["eccentricity"]: - rel["groupName2"] = ( - self.getGroupName(rel["relatedConnection"]) - + "_0DC_" - + self.getGroupName(rel["relatingElement"]) - ) - rel["index"] = len(conn["relatedElements"]) + 1 - rel["unifiedGroupName"] = ( - self.getGroupName(rel["relatedConnection"]) - + "_0DC_%g" % rel["index"] - ) - else: - rel["groupName2"] = self.getGroupName(rel["relatedConnection"]) - self.calculateConstraints(rel) - conn["relatedElements"].append(rel) - # End <-- - - materials = data["db"]["materials"] - profiles = data["db"]["profiles"] - - edgeGroupNames = tuple( - [ - self.getGroupName(el["referenceName"]) - for el in elements - if el["geometryType"] == "line" - ] - ) - faceGroupNames = tuple( - [ - self.getGroupName(el["referenceName"]) - for el in elements - if el["geometryType"] == "surface" - ] - ) - point0DGroupNames = tuple( - [ - self.getGroupName(el["referenceName"]) + "_0D" - for el in connections - if el["geometryType"] == "point" - ] - ) - if includeZeroLength1DSprings: - spring1DGroupNames = tuple( - flatten( - [ - [ - rel["springGroupName"] - for rel in el["connections"] - if rel["springGroupName"] - ] - for el in elements - ] - ) - ) - point1DGroupNames = tuple( - [ - self.getGroupName(el["referenceName"]) + "_0D" - for el in connections - if el["geometryType"] == "line" - ] - ) - - unifiedConnection = False - rigidLinkGroupNames = [] - for conn in connections: - conn["unifiedGroupNames"] = [ - rel["unifiedGroupName"] - for rel in conn["relatedElements"] - if rel["eccentricity"] - ] - # if not conn['appliedCondition'] and len(conn['unifiedGroupNames']) == 1: - # conn['appliedCondition'] = { - # 'dx': True, - # 'dy': True, - # 'dz': True - # } - if len(conn["unifiedGroupNames"]) >= 1: - conn["unifiedGroupNames"].insert(0, self.getGroupName(conn["referenceName"])) - conn["unifiedGroupNames"] = tuple(conn["unifiedGroupNames"]) - unifiedConnection = True - rigidLinkGroupNames.extend( - [ - self.getGroupName(rel["relatingElement"]) - + "_1DR_" - + self.getGroupName(conn["referenceName"]) - for rel in conn["relatedElements"] - if rel["eccentricity"] - ] - ) - rigidLinkGroupNames = tuple(rigidLinkGroupNames) - - # Define file to write command file for code_aster - f = open(self.asterFilename, "w") - - f.write("# Command file generated by IfcOpenShell/ifc2ca scripts\n") - f.write("\n") - - f.write("# Linear Static Analysis With Self-Weight\n") - - f.write( - """ -# STEP: INITIALIZE STUDY -DEBUT( - PAR_LOT = 'NON' -) -""" - ) - - f.write( - """ -# STEP: READ MED FILE -mesh = LIRE_MAILLAGE( - FORMAT = 'MED', - UNITE = 20 -) -""" - ) - - f.write( - """ -# STEP: DEFINE MODEL -model = AFFE_MODELE( - MAILLAGE = mesh, - AFFE = ( - _F( - TOUT = 'OUI', - PHENOMENE = 'MECANIQUE', - MODELISATION = '3D' - ),""" - ) - - if faceGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - PHENOMENE = 'MECANIQUE', - MODELISATION = 'DKT' - ),""" - - context = {"groupNames": faceGroupNames} - - f.write(template.format(**context)) - - if edgeGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - PHENOMENE = 'MECANIQUE', - MODELISATION = 'POU_D_E' - ),""" - - context = {"groupNames": edgeGroupNames} - - f.write(template.format(**context)) - - if point0DGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - PHENOMENE = 'MECANIQUE', - MODELISATION = 'DIS_TR' - ),""" - - context = { - "groupNames": tuple( - flatten( - [ - point0DGroupNames, - spring1DGroupNames if includeZeroLength1DSprings else [], - ] - ) - ) - } - - f.write(template.format(**context)) - - if point1DGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - PHENOMENE = 'MECANIQUE', - MODELISATION = 'DIS_TR' - ),""" - - context = {"groupNames": point1DGroupNames} - - f.write(template.format(**context)) - - if rigidLinkGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - PHENOMENE = 'MECANIQUE', - MODELISATION = 'POU_D_E' - ),""" - - context = {"groupNames": rigidLinkGroupNames} - - f.write(template.format(**context)) - - f.write( - """ - ) -)\n -""" - ) - - f.write("# STEP: DEFINE MATERIALS") - - for i, material in enumerate(materials): - template = """ -{matNameID} = DEFI_MATERIAU( - ELAS = _F( - E = {youngModulus}, - NU = {poissonRatio}, - RHO = {massDensity} - ) -) -""" - if "poissonRatio" in material["mechProps"]: - poissonRatio = material["mechProps"]["poissonRatio"] - else: - if "shearModulus" in material["mechProps"]: - poissonRatio = ( - material["mechProps"]["youngModulus"] - / 2.0 - / material["mechProps"]["shearModulus"] - ) - 1 - else: - poissonRatio = 0.0 - - context = { - "matNameID": "mat" + "_%s" % i, - "youngModulus": float(material["mechProps"]["youngModulus"]), - "poissonRatio": float(poissonRatio), - "massDensity": float(material["commonProps"]["massDensity"]), - } - - f.write(template.format(**context)) - - f.write( - """ -material = AFFE_MATERIAU( - MAILLAGE = mesh, - AFFE = (""" - ) - - for i, material in enumerate(materials): - template = """ - _F( - GROUP_MA = {groupNames}, - MATER = {matNameID}, - ),""" - - context = { - "groupNames": tuple( - [self.getGroupName(rel) for rel in material["relatedElements"]] - ), - "matNameID": "mat" + "_%s" % i, - } - - f.write(template.format(**context)) - - if rigidLinkGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - MATER = {matNameID}, - ),""" - - context = {"groupNames": rigidLinkGroupNames, "matNameID": "mat_0"} - - f.write(template.format(**context)) - - f.write( - """ - ) -) -""" - ) - - f.write( - """ -# STEP: DEFINE ELEMENTS -element = AFFE_CARA_ELEM( - MODELE = model, - POUTRE = (""" - ) - - for profile in profiles: - if ( - profile["profileShape"] == "rectangular" - and profile["profileType"] == "AREA" - ): - template = """ - _F( - GROUP_MA = {groupNames}, - SECTION = 'RECTANGLE', - CARA = ('HY', 'HZ'), - VALE = {profileDimensions} - ),""" - - context = { - "groupNames": tuple( - [self.getGroupName(rel) for rel in profile["relatedElements"]] - ), - "profileDimensions": (profile["xDim"], profile["yDim"]), - } - - f.write(template.format(**context)) - - elif ( - profile["profileShape"] == "iSymmetrical" - and profile["profileType"] == "AREA" - ): - template = """ - _F( - GROUP_MA = {groupNames}, - SECTION = 'GENERALE', - CARA = ('A', 'IY', 'IZ', 'JX'), - VALE = {profileProperties} - ),""" - - context = { - "groupNames": tuple( - [self.getGroupName(rel) for rel in profile["relatedElements"]] - ), - "profileProperties": ( - profile["mechProps"]["crossSectionArea"], - profile["mechProps"]["momentOfInertiaY"], - profile["mechProps"]["momentOfInertiaZ"], - profile["mechProps"]["torsionalConstantX"], - ), - } - - f.write(template.format(**context)) - - if rigidLinkGroupNames: - template = """ - _F( - GROUP_MA = {groupNames}, - SECTION = 'RECTANGLE', - CARA = ('HY', 'HZ'), - VALE = {profileDimensions} - ),""" - - context = {"groupNames": rigidLinkGroupNames, "profileDimensions": (1, 1)} - - f.write(template.format(**context)) - - f.write( - """ - ), - COQUE = (""" - ) - - for el in [el for el in elements if el["geometryType"] == "surface"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - EPAIS = {thickness}, - VECTEUR = {localAxisX} - ),""" - - context = { - "groupName": self.getGroupName(el["referenceName"]), - "thickness": el["thickness"], - "localAxisX": tuple(el["orientation"][0]), - } - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - f.write( - """ - DISCRET = (""" - ) - - for conn in [conn for conn in connections if conn["geometryType"] == "point"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'K_TR_D_N', - VALE = {stiffnesses}, - REPERE = 'LOCAL' - ),""" - - context = { - "groupName": self.getGroupName(conn["referenceName"]) + "_0D", - "stiffnesses": conn["stiffnesses"], - } - - f.write(template.format(**context)) - - if includeZeroLength1DSprings: - for rel in conn["relatedElements"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'K_TR_D_L', - VALE = {stiffnesses}, - REPERE = 'LOCAL' - ),""" - - context = { - "groupName": rel["springGroupName"], - "stiffnesses": rel["stiffnesses"], - } - - f.write(template.format(**context)) - - for conn in [conn for conn in connections if conn["geometryType"] == "line"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'K_TR_D_N', - VALE = {stiffnesses}, - REPERE = 'LOCAL' - ),""" - - context = { - "groupName": self.getGroupName(conn["referenceName"]) + "_0D", - "stiffnesses": conn["stiffnesses"], - } - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - - f.write( - """ - ORIENTATION = (""" - ) - - for el in [el for el in elements if el["geometryType"] == "line"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'VECT_Y', - VALE = {localAxisY} - ),""" - - context = { - "groupName": self.getGroupName(el["referenceName"]), - "localAxisY": tuple(el["orientation"][1]), - } - - f.write(template.format(**context)) - - for conn in [conn for conn in connections if conn["geometryType"] == "point"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'VECT_X_Y', - VALE = {localAxesXY} - ),""" - - context = { - "groupName": self.getGroupName(conn["referenceName"]) + "_0D", - "localAxesXY": tuple(conn["orientation"][0] + conn["orientation"][1]), - } - - f.write(template.format(**context)) - - if includeZeroLength1DSprings: - for rel in conn["relatedElements"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'VECT_X_Y', - VALE = {localAxesXY} - ),""" - - context = { - "groupName": rel["springGroupName"], - "localAxesXY": tuple( - rel["orientation"][0] + rel["orientation"][1] - ), - } - - f.write(template.format(**context)) - - for conn in [conn for conn in connections if conn["geometryType"] == "line"]: - - template = """ - _F( - GROUP_MA = '{groupName}', - CARA = 'VECT_X_Y', - VALE = {localAxesXY} - ),""" - - context = { - "groupName": self.getGroupName(conn["referenceName"]) + "_0D", - "localAxesXY": tuple(conn["orientation"][0] + conn["orientation"][1]), - } - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - - f.write( - """ -)\n -""" - ) - - f.write("# STEP: DEFINE SUPPORTS AND CONSTRAINTS") - - f.write( - """ -liaisons = AFFE_CHAR_MECA( - MODELE = model, - LIAISON_DDL = (""" - ) - - for conn in [conn for conn in connections if conn["geometryType"] == "point"]: - if conn["appliedCondition"]: - for i in range(len(conn["liaisons"]["coeffs"])): - template = """ - _F( - GROUP_NO = {groupNames}, - DDL = {dofs}, - COEF_MULT = {coeffs}, - COEF_IMPO = 0.0 - ),""" - - context = { - "groupNames": conn["liaisons"]["groupNames"], - "dofs": conn["liaisons"]["dofs"][i], - "coeffs": conn["liaisons"]["coeffs"][i], - } - - f.write(template.format(**context)) - - for rel in conn["relatedElements"]: - for i in range(len(rel["liaisons"]["coeffs"])): - template = """ - _F( - GROUP_NO = {groupNames}, - DDL = {dofs}, - COEF_MULT = {coeffs}, - COEF_IMPO = 0.0 - ),""" - - context = { - "groupNames": rel["liaisons"]["groupNames"], - "dofs": rel["liaisons"]["dofs"][i], - "coeffs": rel["liaisons"]["coeffs"][i], - } - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - - f.write( - """ - LIAISON_GROUP = (""" - ) - - for conn in [conn for conn in connections if conn["geometryType"] == "line"]: - if conn["appliedCondition"]: - for i in range(len(conn["liaisons"]["coeffs"])): - template = """ - _F( - GROUP_NO_1 = {groupName_1}, - GROUP_NO_2 = {groupName_1}, - DDL_1 = {dofs}, - DDL_2 = {dofs}, - COEF_MULT_1 = {coeffs}, - COEF_MULT_2 = (0.0, 0.0, 0.0), - COEF_IMPO = 0.0 - ),""" - - context = { - "groupName_1": tuple([conn["liaisons"]["groupNames"][0]]), - "dofs": conn["liaisons"]["dofs"][i], - "coeffs": conn["liaisons"]["coeffs"][i], - } - - f.write(template.format(**context)) - - for rel in conn["relatedElements"]: - for i in range(len(rel["liaisons"]["coeffs"])): - template = """ - _F( - GROUP_NO_1 = {groupName_1}, - GROUP_NO_2 = {groupName_2}, - DDL_1 = {dofs}, - DDL_2 = {dofs}, - COEF_MULT_1 = {coeffs_1}, - COEF_MULT_2 = {coeffs_2}, - COEF_IMPO = 0.0 - ),""" - - context = { - "groupName_1": tuple([rel["liaisons"]["groupNames"][0]]), - "groupName_2": tuple([rel["liaisons"]["groupNames"][3]]), - "dofs": tuple(list(rel["liaisons"]["dofs"][i])[:3]), - "coeffs_1": tuple(list(rel["liaisons"]["coeffs"][i])[:3]), - "coeffs_2": tuple(list(rel["liaisons"]["coeffs"][i])[3:]), - } - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - - if unifiedConnection: - f.write( - """ - LIAISON_UNIF = (""" - ) - - for conn in [ - conn for conn in connections if len(conn["unifiedGroupNames"]) > 1 - ]: - template = """ - _F( - GROUP_NO = {groupNames}, - DDL = ('DX', 'DY', 'DZ', 'DRX', 'DRY', 'DRZ') - ),""" - - context = {"groupNames": conn["unifiedGroupNames"]} - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - - if rigidLinkGroupNames: - f.write( - """ - LIAISON_SOLIDE = (""" - ) - - for groupName in rigidLinkGroupNames: - template = """ - _F( - GROUP_MA = '{groupName}' - ),""" - - context = {"groupName": groupName} - - f.write(template.format(**context)) - - f.write( - """ - ),""" - ) - - f.write( - """ -) -""" - ) - - template = """ -# STEP: DEFINE LOAD -gravLoad = AFFE_CHAR_MECA( - MODELE = model, - PESANTEUR = _F( - GRAVITE = {AccelOfGravity}, - DIRECTION = (0.0, 0.0, -1.0) - ) -) -""" - context = { - "AccelOfGravity": AccelOfGravity, - } - - f.write(template.format(**context)) - - f.write( - """ -# STEP: RUN ANALYSIS -res_Bld = MECA_STATIQUE( - MODELE = model, - CHAM_MATER = material, - CARA_ELEM = element, - EXCIT = ( - _F( - CHARGE = liaisons - ), - _F( - CHARGE = gravLoad - ) - ) -) -""" - ) - - # f.write( - # ''' - # # STEP: POST-PROCESSING - # res_Bld = CALC_CHAMP( - # reuse = res_Bld, - # RESULTAT = res_Bld, - # # CONTRAINTE = ('SIEF_ELNO', 'SIGM_ELNO', 'EFGE_ELNO',), - # FORCE = ('REAC_NODA', 'FORC_NODA',) - # ) - # ''' - # ) - # - # template = \ - # ''' - # # STEP: MASS EXTRACTION FOR EACH ASSEMBLE - # FaceMass = POST_ELEM( - # TITRE = 'TotMass', - # MODELE = model, - # CARA_ELEM = element, - # CHAM_MATER = material, - # MASS_INER = _F( - # GROUP_MA = {massList}, - # ), - # )\n''' - # - # context = { - # 'massList': massList, - # } - # - # f.write(template.format(**context)) - # - # f.write( - # ''' - # IMPR_TABLE( - # UNITE = 10, - # TABLE = FaceMass, - # SEPARATEUR = ',', - # NOM_PARA = ('LIEU', 'MASSE', 'CDG_X', 'CDG_Y', 'CDG_Z'), - # # FORMAT_R = '1PE15.6', - # ) - # ''' - # ) - # - # template = \ - # ''' - # # STEP: REACTION EXTRACTION AT THE BASE - # Reacs = POST_RELEVE_T( - # ACTION = _F( - # INTITULE = 'sumReac', - # GROUP_NO = {groupNames}, - # RESULTAT = res_Bld, - # NOM_CHAM = 'REAC_NODA', - # RESULTANTE = ('DX','DY','DZ',), - # MOMENT = ('DRX','DRY','DRZ',), - # POINT = (0,0,0,), - # OPERATION = 'EXTRACTION' - # ) - # ) - # ''' - # - # context = { - # 'groupNames': point0DGroupNames, - # } - # - # f.write(template.format(**context)) - # - # f.write( - # ''' - # IMPR_TABLE( - # UNITE = 10, - # TABLE = Reacs, - # SEPARATEUR = ',', - # # NOM_PARA = ('INTITULE', 'RESU', 'NOM_CHAM', 'INST', 'DX','DY','DZ'), - # FORMAT_R = '1PE12.3', - # ) - # ''' - # ) - # - f.write( - """ -# STEP: DEFORMED SHAPE EXTRACTION -IMPR_RESU( - FORMAT = 'MED', - UNITE = 80, - RESU = _F( - RESULTAT = res_Bld, - NOM_CHAM = ('DEPL',), # 'REAC_NODA', 'FORC_NODA', - NOM_CHAM_MED = ('Bld_DISP',), # 'Bld_REAC', 'Bld_FORC' - ) -) -""" - ) - - f.write( - """ -# STEP: CONCLUDE STUDY -FIN() -""" - ) - - f.close() - def calculateConstraints(self, rel): gr1 = rel["groupName1"] gr2 = rel["groupName2"] @@ -914,95 +61,47 @@ FIN() "dry": True, "drz": True, } - if ( - isinstance(rel["appliedCondition"]["dx"], bool) - and rel["appliedCondition"]["dx"] - ): - liaisons["coeffs"].append( - (o[0][0], o[1][0], o[2][0], -o[0][0], -o[1][0], -o[2][0]) - ) + if isinstance(rel["appliedCondition"]["dx"], bool) and rel["appliedCondition"]["dx"]: + liaisons["coeffs"].append((o[0][0], o[1][0], o[2][0], -o[0][0], -o[1][0], -o[2][0])) liaisons["dofs"].append(("DX", "DY", "DZ", "DX", "DY", "DZ")) - elif ( - isinstance(rel["appliedCondition"]["dx"], float) - and rel["appliedCondition"]["dx"] > 0 - ): + elif isinstance(rel["appliedCondition"]["dx"], float) and rel["appliedCondition"]["dx"] > 0: stiffnesses[0] = rel["appliedCondition"]["dx"] - if ( - isinstance(rel["appliedCondition"]["dy"], bool) - and rel["appliedCondition"]["dy"] - ): - liaisons["coeffs"].append( - (o[0][1], o[1][1], o[2][1], -o[0][1], -o[1][1], -o[2][1]) - ) + if isinstance(rel["appliedCondition"]["dy"], bool) and rel["appliedCondition"]["dy"]: + liaisons["coeffs"].append((o[0][1], o[1][1], o[2][1], -o[0][1], -o[1][1], -o[2][1])) liaisons["dofs"].append(("DX", "DY", "DZ", "DX", "DY", "DZ")) - elif ( - isinstance(rel["appliedCondition"]["dy"], float) - and rel["appliedCondition"]["dy"] > 0 - ): + elif isinstance(rel["appliedCondition"]["dy"], float) and rel["appliedCondition"]["dy"] > 0: stiffnesses[1] = rel["appliedCondition"]["dy"] - if ( - isinstance(rel["appliedCondition"]["dz"], bool) - and rel["appliedCondition"]["dz"] - ): - liaisons["coeffs"].append( - (o[0][2], o[1][2], o[2][2], -o[0][2], -o[1][2], -o[2][2]) - ) + if isinstance(rel["appliedCondition"]["dz"], bool) and rel["appliedCondition"]["dz"]: + liaisons["coeffs"].append((o[0][2], o[1][2], o[2][2], -o[0][2], -o[1][2], -o[2][2])) liaisons["dofs"].append(("DX", "DY", "DZ", "DX", "DY", "DZ")) - elif ( - isinstance(rel["appliedCondition"]["dz"], float) - and rel["appliedCondition"]["dz"] > 0 - ): + elif isinstance(rel["appliedCondition"]["dz"], float) and rel["appliedCondition"]["dz"] > 0: stiffnesses[2] = rel["appliedCondition"]["dz"] - if ( - isinstance(rel["appliedCondition"]["drx"], bool) - and rel["appliedCondition"]["drx"] - ): - liaisons["coeffs"].append( - (o[0][0], o[1][0], o[2][0], -o[0][0], -o[1][0], -o[2][0]) - ) + if isinstance(rel["appliedCondition"]["drx"], bool) and rel["appliedCondition"]["drx"]: + liaisons["coeffs"].append((o[0][0], o[1][0], o[2][0], -o[0][0], -o[1][0], -o[2][0])) liaisons["dofs"].append(("DRX", "DRY", "DRZ", "DRX", "DRY", "DRZ")) - elif ( - isinstance(rel["appliedCondition"]["drx"], float) - and rel["appliedCondition"]["drx"] > 0 - ): + elif isinstance(rel["appliedCondition"]["drx"], float) and rel["appliedCondition"]["drx"] > 0: stiffnesses[3] = rel["appliedCondition"]["drx"] - if ( - isinstance(rel["appliedCondition"]["dry"], bool) - and rel["appliedCondition"]["dry"] - ): - liaisons["coeffs"].append( - (o[0][1], o[1][1], o[2][1], -o[0][1], -o[1][1], -o[2][1]) - ) + if isinstance(rel["appliedCondition"]["dry"], bool) and rel["appliedCondition"]["dry"]: + liaisons["coeffs"].append((o[0][1], o[1][1], o[2][1], -o[0][1], -o[1][1], -o[2][1])) liaisons["dofs"].append(("DRX", "DRY", "DRZ", "DRX", "DRY", "DRZ")) - elif ( - isinstance(rel["appliedCondition"]["dry"], float) - and rel["appliedCondition"]["dry"] > 0 - ): + elif isinstance(rel["appliedCondition"]["dry"], float) and rel["appliedCondition"]["dry"] > 0: stiffnesses[4] = rel["appliedCondition"]["dry"] - if ( - isinstance(rel["appliedCondition"]["drz"], bool) - and rel["appliedCondition"]["drz"] - ): - liaisons["coeffs"].append( - (o[0][2], o[1][2], o[2][2], -o[0][2], -o[1][2], -o[2][2]) - ) + if isinstance(rel["appliedCondition"]["drz"], bool) and rel["appliedCondition"]["drz"]: + liaisons["coeffs"].append((o[0][2], o[1][2], o[2][2], -o[0][2], -o[1][2], -o[2][2])) liaisons["dofs"].append(("DRX", "DRY", "DRZ", "DRX", "DRY", "DRZ")) - elif ( - isinstance(rel["appliedCondition"]["drz"], float) - and rel["appliedCondition"]["drz"] > 0 - ): + elif isinstance(rel["appliedCondition"]["drz"], float) and rel["appliedCondition"]["drz"] > 0: stiffnesses[5] = rel["appliedCondition"]["drz"] rel["liaisons"] = liaisons rel["stiffnesses"] = tuple(stiffnesses) def calculateRestraints(self, conn): - group = self.getGroupName(conn["referenceName"]) + group = self.getGroupName(conn["ref_id"]) o = np.array(conn["orientation"]).transpose().tolist() liaisons = {"groupNames": (group, group, group), "coeffs": [], "dofs": []} stiffnesses = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] @@ -1012,96 +111,299 @@ FIN() conn["stiffnesses"] = tuple(stiffnesses) return - if ( - isinstance(conn["appliedCondition"]["dx"], bool) - and conn["appliedCondition"]["dx"] - ): + if isinstance(conn["appliedCondition"]["dx"], bool) and conn["appliedCondition"]["dx"]: liaisons["coeffs"].append((o[0][0], o[1][0], o[2][0])) liaisons["dofs"].append(("DX", "DY", "DZ")) - elif ( - isinstance(conn["appliedCondition"]["dx"], float) - and conn["appliedCondition"]["dx"] > 0 - ): + elif isinstance(conn["appliedCondition"]["dx"], float) and conn["appliedCondition"]["dx"] > 0: stiffnesses[0] = conn["appliedCondition"]["dx"] - if ( - isinstance(conn["appliedCondition"]["dy"], bool) - and conn["appliedCondition"]["dy"] - ): + if isinstance(conn["appliedCondition"]["dy"], bool) and conn["appliedCondition"]["dy"]: liaisons["coeffs"].append((o[0][1], o[1][1], o[2][1])) liaisons["dofs"].append(("DX", "DY", "DZ")) - elif ( - isinstance(conn["appliedCondition"]["dy"], float) - and conn["appliedCondition"]["dy"] > 0 - ): + elif isinstance(conn["appliedCondition"]["dy"], float) and conn["appliedCondition"]["dy"] > 0: stiffnesses[1] = conn["appliedCondition"]["dy"] - if ( - isinstance(conn["appliedCondition"]["dz"], bool) - and conn["appliedCondition"]["dz"] - ): + if isinstance(conn["appliedCondition"]["dz"], bool) and conn["appliedCondition"]["dz"]: liaisons["coeffs"].append((o[0][2], o[1][2], o[2][2])) liaisons["dofs"].append(("DX", "DY", "DZ")) - elif ( - isinstance(conn["appliedCondition"]["dz"], float) - and conn["appliedCondition"]["dz"] > 0 - ): + elif isinstance(conn["appliedCondition"]["dz"], float) and conn["appliedCondition"]["dz"] > 0: stiffnesses[2] = conn["appliedCondition"]["dz"] - if ( - isinstance(conn["appliedCondition"]["drx"], bool) - and conn["appliedCondition"]["drx"] - ): + if isinstance(conn["appliedCondition"]["drx"], bool) and conn["appliedCondition"]["drx"]: liaisons["coeffs"].append((o[0][0], o[1][0], o[2][0])) liaisons["dofs"].append(("DRX", "DRY", "DRZ")) - elif ( - isinstance(conn["appliedCondition"]["drx"], float) - and conn["appliedCondition"]["drx"] > 0 - ): + elif isinstance(conn["appliedCondition"]["drx"], float) and conn["appliedCondition"]["drx"] > 0: stiffnesses[3] = conn["appliedCondition"]["drx"] - if ( - isinstance(conn["appliedCondition"]["dry"], bool) - and conn["appliedCondition"]["dry"] - ): + if isinstance(conn["appliedCondition"]["dry"], bool) and conn["appliedCondition"]["dry"]: liaisons["coeffs"].append((o[0][1], o[1][1], o[2][1])) liaisons["dofs"].append(("DRX", "DRY", "DRZ")) - elif ( - isinstance(conn["appliedCondition"]["dry"], float) - and conn["appliedCondition"]["dry"] > 0 - ): + elif isinstance(conn["appliedCondition"]["dry"], float) and conn["appliedCondition"]["dry"] > 0: stiffnesses[4] = conn["appliedCondition"]["dry"] - if ( - isinstance(conn["appliedCondition"]["drz"], bool) - and conn["appliedCondition"]["drz"] - ): + if isinstance(conn["appliedCondition"]["drz"], bool) and conn["appliedCondition"]["drz"]: liaisons["coeffs"].append((o[0][2], o[1][2], o[2][2])) liaisons["dofs"].append(("DRX", "DRY", "DRZ")) - elif ( - isinstance(conn["appliedCondition"]["drz"], float) - and conn["appliedCondition"]["drz"] > 0 - ): + elif isinstance(conn["appliedCondition"]["drz"], float) and conn["appliedCondition"]["drz"] > 0: stiffnesses[5] = conn["appliedCondition"]["drz"] conn["liaisons"] = liaisons conn["stiffnesses"] = tuple(stiffnesses) + def create_comm(self, path, cases): + self.comm_path = Path(path) + self.cases = cases -if __name__ == "__main__": - fileNames = [ - "cantilever_01", - "portal_01", - "grid_of_beams", - "slab_01", - "structure_01", - ] - files = fileNames + data = self.data - for fileName in files: - BASE_PATH = Path( - "/home/jesusbill/Dev-Projects/github.com/IfcOpenShell/analysis-models/models/" + elements = data["elements"] + connections = data["connections"] + # --> Delete this reference data and repopulate it with the objects + # while going through elements + for conn in connections: + conn["related_elements"] = [] + self.calculateRestraints(conn) + for el in elements: + for rel in el["connections"]: + conn = [c for c in connections if c["ref_id"] == rel["related_connection"]][0] + rel["conn_string"] = None + if conn["geometry_type"] == "Vertex": + rel["conn_string"] = "_0DC_" + rel["springGroupName"] = ( + self.getGroupName(rel["relating_element"]) + + "_1DS_" + + self.getGroupName(rel["related_connection"]) + ) + if conn["geometry_type"] == "Edge": + rel["conn_string"] = "_1DC_" + rel["springGroupName"] = None + if conn["geometry_type"] == "Face": + rel["conn_string"] = "_2DC_" + rel["springGroupName"] = None + + rel["groupName1"] = ( + self.getGroupName(rel["relating_element"]) + + rel["conn_string"] + + self.getGroupName(rel["related_connection"]) + ) + if rel["eccentricity"]: + rel["groupName2"] = ( + self.getGroupName(rel["related_connection"]) + + "_0DC_" + + self.getGroupName(rel["relating_element"]) + ) + rel["index"] = len(conn["related_elements"]) + 1 + rel["unifiedGroupName"] = self.getGroupName(rel["related_connection"]) + "_0DC_%g" % rel["index"] + else: + rel["groupName2"] = self.getGroupName(rel["related_connection"]) + self.calculateConstraints(rel) + conn["related_elements"].append(rel) + # End <-- + + materials = data["db"]["materials"] + for _, material in materials.items(): + material["groupNames"] = tuple([self.getGroupName(rel) for rel in material["related_elements"]]) + + profiles = data["db"]["profiles"] + for _, profile in profiles.items(): + profile["groupNames"] = tuple([self.getGroupName(rel) for rel in profile["related_elements"]]) + + edgeGroupNames = tuple([self.getGroupName(el["ref_id"]) for el in elements if el["geometry_type"] == "Edge"]) + faceGroupNames = tuple([self.getGroupName(el["ref_id"]) for el in elements if el["geometry_type"] == "Face"]) + point0DGroupNames = tuple( + [self.getGroupName(el["ref_id"]) + "_0D" for el in connections if el["geometry_type"] == "Vertex"] ) - DATAFILENAME = BASE_PATH / fileName / f"{fileName}.json" - ASTERFILENAME = BASE_PATH / fileName / f"{fileName}.comm" - COMMANDFILE(DATAFILENAME, ASTERFILENAME) + if includeZeroLength1DSprings: + spring1DGroupNames = tuple( + flatten( + [[rel["springGroupName"] for rel in el["connections"] if rel["springGroupName"]] for el in elements] + ) + ) + point1DGroupNames = tuple( + [self.getGroupName(el["ref_id"]) + "_0D" for el in connections if el["geometry_type"] == "Edge"] + ) + + rigidLinkGroupNames = [] + for conn in connections: + conn["unifiedGroupNames"] = [ + rel["unifiedGroupName"] for rel in conn["related_elements"] if rel["eccentricity"] + ] + # if not conn['appliedCondition'] and len(conn['unifiedGroupNames']) == 1: + # conn['appliedCondition'] = { + # 'dx': True, + # 'dy': True, + # 'dz': True + # } + if len(conn["unifiedGroupNames"]): + conn["unifiedGroupNames"].insert(0, self.getGroupName(conn["ref_id"])) + conn["unifiedGroupNames"] = tuple(conn["unifiedGroupNames"]) + rigidLinkGroupNames.extend( + [ + self.getGroupName(rel["relating_element"]) + "_1DR_" + self.getGroupName(conn["ref_id"]) + for rel in conn["related_elements"] + if rel["eccentricity"] + ] + ) + rigidLinkGroupNames = tuple(rigidLinkGroupNames) + + # Start of Writng Command File# Define file to write command file for code_aster + with self.comm_path.open("w") as f: + start_process = self.env.get_template("codeaster/start_process.py") + f.write(start_process.render()) + + read_mesh = self.env.get_template("codeaster/read_mesh.py") + f.write(read_mesh.render()) + + define_model = self.env.get_template("codeaster/define_model.py") + f.write( + define_model.render( + faceGroupNames=faceGroupNames, + edgeGroupNames=edgeGroupNames, + point0DGroupNames=point0DGroupNames, + point0DGroupNamesPlus=point0DGroupNames + + (spring1DGroupNames if includeZeroLength1DSprings else tuple()), + point1DGroupNames=point1DGroupNames, + rigidLinkGroupNames=rigidLinkGroupNames, + ) + ) + + define_materials = self.env.get_template("codeaster/define_materials.py") + f.write( + define_materials.render( + enumerate=enumerate, + materials=materials, + rigidLinkGroupNames=rigidLinkGroupNames, + ) + ) + + define_elements = self.env.get_template("codeaster/define_elements.py") + f.write( + define_elements.render( + len=len, + tuple=tuple, + getGroupName=self.getGroupName, + includeZeroLength1DSprings=includeZeroLength1DSprings, + profiles=profiles, + rigidLinkGroupNames=rigidLinkGroupNames, + beamElements=[el for el in elements if el["geometry_type"] == "Edge"], + shellElements=[el for el in elements if el["geometry_type"] == "Face"], + vertexConnections=[conn for conn in connections if conn["geometry_type"] == "Vertex"], + edgeConnections=[conn for conn in connections if conn["geometry_type"] == "Edge"], + # faceConnections=[conn for conn in connections if conn["geometry_type"] == "Face"], + ) + ) + + define_connections = self.env.get_template("codeaster/define_connections.py") + f.write( + define_connections.render( + len=len, + range=range, + tuple=tuple, + rigidLinkGroupNames=rigidLinkGroupNames, + vertexConnections=[conn for conn in connections if conn["geometry_type"] == "Vertex"], + edgeConnections=[conn for conn in connections if conn["geometry_type"] == "Edge"], + unifiedConnections=[conn for conn in connections if len(conn["unifiedGroupNames"]) > 1], + ) + ) + + for case_instant in self.cases: + if case_instant == "LC": + define_loads = self.env.get_template("codeaster/define_loads.py") + f.write( + define_loads.render( + tuple=tuple, + analysis_time=f"analysis_time_{case_instant}", + load=f"load_{case_instant}", + load_key=f"loads{case_instant}", + start=1.0, + end=float(len(self.data["load_cases"])), + steps=len(self.data["load_cases"]) - 1, + getGroupName=self.getGroupName, + time=tuple([float(t) for t in (range(1, len(self.data["load_cases"]) + 1))]), + vertexLoadElements=[ + item + for item in connections + if item["geometry_type"] == "Vertex" and item["loads"] is not None + ], + edgeLoadElements=[ + item + for item in elements + connections + if item["geometry_type"] == "Edge" and item["loads"] is not None + ], + faceLoadElements=[ + item + for item in elements + connections + if item["geometry_type"] == "Face" and item["loads"] is not None + ], + ) + ) + + elif case_instant == "COMB": + define_loads = self.env.get_template("codeaster/define_loads.py") + f.write( + define_loads.render( + tuple=tuple, + analysis_time=f"analysis_time_{case_instant}", + load=f"load_{case_instant}", + load_key=f"loads{case_instant}", + start=1.0, + end=float(len(self.data["load_combinations"])), + steps=len(self.data["load_combinations"]) - 1, + getGroupName=self.getGroupName, + time=tuple([float(t) for t in (range(1, len(self.data["load_combinations"]) + 1))]), + vertexLoadElements=[ + item + for item in connections + if item["geometry_type"] == "Vertex" and item["loads"] is not None + ], + edgeLoadElements=[ + item + for item in elements + connections + if item["geometry_type"] == "Edge" and item["loads"] is not None + ], + faceLoadElements=[ + item + for item in elements + connections + if item["geometry_type"] == "Face" and item["loads"] is not None + ], + ) + ) + + for case_instant in self.cases: + run_analysis = self.env.get_template("codeaster/run_analysis.py") + if case_instant == "LC": + f.write( + run_analysis.render( + analysis_time="analysis_time_LC", + load=f"load_{case_instant}", + res_Bld=f"res_Bld_{case_instant}", + ) + ) + elif case_instant == "COMB": + f.write( + run_analysis.render( + analysis_time="analysis_time_COMB", + load=f"load_{case_instant}", + res_Bld=f"res_Bld_{case_instant}", + ) + ) + + for case_instant in self.cases: + export_results = self.env.get_template("codeaster/export_results.py") + if case_instant == "LC": + f.write( + export_results.render( + unit_number=80, + res_Bld=f"res_Bld_{case_instant}", + ) + ) + elif case_instant == "COMB": + f.write( + export_results.render( + unit_number=81, + res_Bld=f"res_Bld_{case_instant}", + ) + ) + + finish_process = self.env.get_template("codeaster/finish_process.py") + f.write(finish_process.render()) diff --git a/src/ifc2ca/templates/codeaster/define_connections.py b/src/ifc2ca/templates/codeaster/define_connections.py new file mode 100644 index 0000000000..6636c44688 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/define_connections.py @@ -0,0 +1,80 @@ +# STEP: DEFINE SUPPORTS AND CONSTRAINTS +connection = AFFE_CHAR_MECA( + MODELE = model, + {%- if vertexConnections %} + LIAISON_DDL = ( + {%- for conn in vertexConnections %} + {%- if conn.appliedCondition %} + {%- for i in range(len(conn.liaisons.coeffs)) %} + _F( + GROUP_NO = {{ conn.liaisons.groupNames }}, + DDL = {{ conn.liaisons.dofs[i] }}, + COEF_MULT = {{ conn.liaisons.coeffs[i] }}, + COEF_IMPO = 0.0 + ), + {%- endfor %} + {%- endif %} + {%- for rel in conn.related_elements %} + {%- for i in range(len(rel.liaisons.coeffs)) %} + _F( + GROUP_NO = {{ rel.liaisons.groupNames }}, + DDL = {{ rel.liaisons.dofs[i] }}, + COEF_MULT = {{ rel.liaisons.coeffs[i] }}, + COEF_IMPO = 0.0 + ), + {%- endfor %} + {%- endfor %} + {%- endfor %} + ), + {%- endif %} + {%- if edgeConnections %} + LIAISON_GROUP = ( + {%- for conn in edgeConnections %} + {%- if conn.appliedCondition %} + {%- for i in range(len(conn.liaisons.coeffs)) %} + _F( + GROUP_NO_1 = {{ tuple([conn.liaisons.groupNames[0]]) }}, + GROUP_NO_2 = {{ tuple([conn.liaisons.groupNames[0]]) }}, + DDL_1 = {{ conn.liaisons.dofs[i] }}, + DDL_2 = {{ conn.liaisons.dofs[i] }}, + COEF_MULT_1 = {{ conn.liaisons.coeffs[i] }}, + COEF_MULT_2 = (0.0, 0.0, 0.0), + COEF_IMPO = 0.0 + ), + {%- endfor %} + {%- endif %} + {%- for rel in conn.related_elements %} + {%- for i in range(len(rel.liaisons.coeffs)) %} + _F( + GROUP_NO_1 = {{ tuple([rel.liaisons.groupNames[0]]) }}, + GROUP_NO_2 = {{ tuple([rel.liaisons.groupNames[3]]) }}, + DDL_1 = {{ tuple(rel.liaisons.dofs[i][:3]) }}, + DDL_2 = {{ tuple(rel.liaisons.dofs[i][:3]) }}, + COEF_MULT_1 = {{ tuple(rel.liaisons.coeffs[i][:3]) }}, + COEF_MULT_2 = {{ tuple(rel.liaisons.coeffs[i][3:]) }}, + COEF_IMPO = 0.0 + ), + {%- endfor %} + {%- endfor %} + {%- endfor %} + ), + {%- endif %} + {%- if unifiedConnections %} + LIAISON_UNIF = ( + {%- for conn in unifiedConnections %} + _F( + GROUP_NO = {{ conn.unifiedGroupNames }}, + DDL = ('DX', 'DY', 'DZ', 'DRX', 'DRY', 'DRZ') + ), + {%- endfor %} + ), + {%- endif %} + {%- if rigidLinkGroupNames %} + LIAISON_SOLIDE = ( + {%- for groupName in rigidLinkGroupNames %} + _F(GROUP_MA = {{ tuple([groupName]) }}), + {%- endfor %} + ), + {%- endif %} +) +{{ "\n" }} diff --git a/src/ifc2ca/templates/codeaster/define_elements.py b/src/ifc2ca/templates/codeaster/define_elements.py new file mode 100644 index 0000000000..3fb3f1e0fa --- /dev/null +++ b/src/ifc2ca/templates/codeaster/define_elements.py @@ -0,0 +1,115 @@ +# STEP: DEFINE ELEMENTS +element = AFFE_CARA_ELEM( + MODELE = model, + POUTRE = ( + {%- for _, profile in profiles.items() %} + {%- if profile.properties %} + _F( + GROUP_MA = {{ profile.groupNames }}, + SECTION = 'GENERALE', + CARA = ('A', 'IY', 'IZ', 'JX'), + VALE = ({{ profile.properties.CrossSectionArea }}, {{ profile.properties.MomentOfInertiaY }}, {{ profile.properties.MomentOfInertiaZ }}, {{ profile.properties.TorsionalConstantX }}) + ), + {%- elif profile.type == "IfcRectangleProfileDef" and profile.ProfileType == "AREA" %} + _F( + GROUP_MA = {{ profile.groupNames }}, + SECTION = 'RECTANGLE', + CARA = ('HY', 'HZ'), + VALE = ({{ profile.XDim }}, {{ profile.YDim }}) + ), + {%- elif profile.type == "IfcRectangleHollowProfileDef" and profile.ProfileType == "AREA" %} + _F( + GROUP_MA = {{ profile.groupNames }}, + SECTION = 'RECTANGLE', + CARA = ('HY', 'HZ', 'EPY', 'EPZ'), + VALE = ({{ profile.XDim }}, {{ profile.YDim }}, {{ profile.WallThickness }}, {{ profile.WallThickness }}) + ), + {%- else %} + _F( + GROUP_MA = {{ profile.groupNames }}, + SECTION = 'GENERALE', + CARA = ('A', 'IY', 'IZ', 'JX'), + VALE = ({{ profile.properties.CrossSectionArea }}, {{ profile.properties.MomentOfInertiaY }}, {{ profile.properties.MomentOfInertiaZ }}, {{ profile.properties.TorsionalConstantX }}) + ), + {%- endif %} + {%- endfor %} + {%- if rigidLinkGroupNames %} + _F( + GROUP_MA = {{ rigidLinkGroupNames }}, + SECTION = 'RECTANGLE', + CARA = ('HY', 'HZ'), + VALE = (1.0, 1.0) + ), + {%- endif %} + ), + COQUE = ( + {%- for el in shellElements %} + _F( + GROUP_MA = {{ tuple([getGroupName(el.ref_id)]) }}, + EPAIS = {{ el.Thickness }}, + VECTEUR = {{ tuple(el.orientation[0]) }} + ), + {%- endfor %} + ), + DISCRET = ( + {%- for conn in vertexConnections %} + _F( + GROUP_MA = {{ tuple([getGroupName(conn.ref_id) + "_0D"]) }}, + CARA = 'K_TR_D_N', + VALE = {{ conn.stiffnesses }}, + REPERE = 'LOCAL' + ), + {%- if includeZeroLength1DSprings %} + {%- for rel in conn.related_elements %} + _F( + GROUP_MA = {{ tuple([rel.springGroupName]) }}, + CARA = 'K_TR_D_L', + VALE = {{ rel.stiffnesses }}, + REPERE = 'LOCAL' + ), + {%- endfor %} + {%- endif %} + {%- endfor %} + {%- for conn in edgeConnections %} + _F( + GROUP_MA = {{ tuple([getGroupName(conn.ref_id) + "_0D"]) }}, + CARA = 'K_TR_D_N', + VALE = {{ conn.stiffnesses }}, + REPERE = 'LOCAL' + ), + {%- endfor %} + ), + ORIENTATION = ( + {%- for el in beamElements %} + _F( + GROUP_MA = {{ tuple([getGroupName(el.ref_id)]) }}, + CARA = 'VECT_Y', + VALE = {{ tuple(el.orientation[1]) }} + ), + {%- endfor %} + {%- for conn in vertexConnections %} + _F( + GROUP_MA = {{ tuple([getGroupName(conn.ref_id) + "_0D"]) }}, + CARA = 'VECT_X_Y', + VALE = {{ tuple(conn.orientation[0] + conn.orientation[1]) }} + ), + {%- if includeZeroLength1DSprings %} + {%- for rel in conn.related_elements %} + _F( + GROUP_MA = {{ tuple([rel.springGroupName]) }}, + CARA = 'VECT_X_Y', + VALE = {{ tuple(rel.orientation[0] + rel.orientation[1]) }}, + ), + {%- endfor %} + {%- endif %} + {%- endfor %} + {%- for conn in edgeConnections %} + _F( + GROUP_MA = {{ tuple([getGroupName(conn.ref_id) + "_0D"]) }}, + CARA = 'VECT_X_Y', + VALE = {{ tuple(conn.orientation[0] + conn.orientation[1]) }} + ), + {%- endfor %} + ), +) +{{ "\n" }} diff --git a/src/ifc2ca/templates/codeaster/define_loads.py b/src/ifc2ca/templates/codeaster/define_loads.py new file mode 100644 index 0000000000..c65bb625b4 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/define_loads.py @@ -0,0 +1,44 @@ +# STEP: DEFINE TIME +{{ analysis_time }} = DEFI_LIST_REEL( + DEBUT = {{ start }}, + INTERVALLE = _F( + JUSQU_A = {{ end }}, + NOMBRE = {{ steps }} + ) +) + +# STEP: DEFINE LOADS +{{ load }} = AFFE_CHAR_MECA_F( + MODELE = model, + FORCE_NODALE = ( + {%- for el in vertexLoadElements %} + _F( + GROUP_NO = {{ tuple([getGroupName(el.ref_id)]) }}, + {%- for key, load in el.loads[load_key].items() %} + {{ key }} = DEFI_FONCTION(NOM_PARA='INST', ABSCISSE={{ time }}, ORDONNEE={{ tuple(load) }}), + {%- endfor %} + ), + {%- endfor %} + ), + FORCE_POUTRE = ( + {%- for el in edgeLoadElements %} + _F( + GROUP_MA = {{ tuple([getGroupName(el.ref_id)]) }}, + {%- for key, load in el.loads[load_key].items() %} + {{ key }} = DEFI_FONCTION(NOM_PARA='INST', ABSCISSE={{ time }}, ORDONNEE={{ tuple(load) }}), + {%- endfor %} + ), + {%- endfor %} + ), + FORCE_COQUE = ( + {%- for el in faceLoadElements %} + _F( + GROUP_MA = {{ tuple([getGroupName(el.ref_id)]) }}, + {%- for key, load in el.loads[load_key].items() %} + {{ key }} = DEFI_FONCTION(NOM_PARA='INST', ABSCISSE={{ time }}, ORDONNEE={{ tuple(load) }}), + {%- endfor %} + ), + {%- endfor %} + ), +) +{{ "\n" }} diff --git a/src/ifc2ca/templates/codeaster/define_materials.py b/src/ifc2ca/templates/codeaster/define_materials.py new file mode 100644 index 0000000000..b4d5b995ae --- /dev/null +++ b/src/ifc2ca/templates/codeaster/define_materials.py @@ -0,0 +1,28 @@ +# STEP: DEFINE MATERIALS +{%- for i, (_, material) in enumerate(materials.items()) %} +{{ "mat" + "_%s" % i }} = DEFI_MATERIAU( + ELAS = _F( + E = {{ material.properties.YoungModulus }}, + NU = {{ material.properties.PoissonRatio }}, + RHO = {{ material.properties.MassDensity }} + ) +) +{% endfor %} +material = AFFE_MATERIAU( + MAILLAGE = mesh, + AFFE = ( + {%- for i, (_, material) in enumerate(materials.items()) %} + _F( + GROUP_MA = {{ material.groupNames }}, + MATER = {{ "mat" + "_%s" % i }}, + ), + {%- endfor %} + {%- if rigidLinkGroupNames %} + _F( + GROUP_MA = {{ rigidLinkGroupNames }}, + MATER = {{ "mat_0" }}, + ), + {%- endif %} + ) +) +{{ "\n" }} diff --git a/src/ifc2ca/templates/codeaster/define_model.py b/src/ifc2ca/templates/codeaster/define_model.py new file mode 100644 index 0000000000..6aea3eec81 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/define_model.py @@ -0,0 +1,47 @@ +# STEP: DEFINE MODEL +model = AFFE_MODELE( + MAILLAGE = mesh, + AFFE = ( + _F( + TOUT = 'OUI', + PHENOMENE = 'MECANIQUE', + MODELISATION = '3D' + ), + {%- if faceGroupNames %} + _F( + GROUP_MA = {{ faceGroupNames }}, + PHENOMENE = 'MECANIQUE', + MODELISATION = 'DKT' + ), + {%- endif %} + {%- if edgeGroupNames %} + _F( + GROUP_MA = {{ edgeGroupNames }}, + PHENOMENE = 'MECANIQUE', + MODELISATION = 'POU_D_E' + ), + {%- endif %} + {%- if point0DGroupNames %} + _F( + GROUP_MA = {{ point0DGroupNamesPlus }}, + PHENOMENE = 'MECANIQUE', + MODELISATION = 'DIS_TR' + ), + {%- endif %} + {%- if point1DGroupNames %} + _F( + GROUP_MA = {{ point1DGroupNames }}, + PHENOMENE = 'MECANIQUE', + MODELISATION = 'DIS_TR' + ), + {%- endif %} + {%- if rigidLinkGroupNames %} + _F( + GROUP_MA = {{ rigidLinkGroupNames }}, + PHENOMENE = 'MECANIQUE', + MODELISATION = 'POU_D_E' + ), + {%- endif %} + ) +) +{{ "\n" }} diff --git a/src/ifc2ca/templates/codeaster/export b/src/ifc2ca/templates/codeaster/export new file mode 100644 index 0000000000..2f7daaf7a2 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/export @@ -0,0 +1,13 @@ +P actions make_etude +P memory_limit {{ allocated_memory }} +P time_limit {{ time_limit }} +P version stable +F comm {{ model_name }}_{{ run_label }}.comm D 1 +F libr {{ model_name }}.med D 20 +F mess {{ model_name }}_{{ run_label }}.mess R 6 +{%- if "LC" in cases %} +F rmed {{ model_name + "_LC" }}.rmed R 80 +{%- endif %} +{%- if "COMB" in cases %} +F rmed {{ model_name + "_COMB" }}.rmed R 81 +{%- endif %} diff --git a/src/ifc2ca/templates/codeaster/export_results.py b/src/ifc2ca/templates/codeaster/export_results.py new file mode 100644 index 0000000000..c5c5265538 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/export_results.py @@ -0,0 +1,11 @@ +# STEP: RESULT EXTRACTION +IMPR_RESU( + FORMAT="MED", + UNITE={{unit_number}}, + RESU=_F( + RESULTAT={{res_Bld}}, + NOM_CHAM=("DEPL", "EFGE_NOEU"), + NOM_CHAM_MED=("MODEL_DISP", "ELEMENT_FORCE"), + ), +) +{{"\n"}} diff --git a/src/ifc2ca/templates/codeaster/finish_process.py b/src/ifc2ca/templates/codeaster/finish_process.py new file mode 100644 index 0000000000..69009c67ca --- /dev/null +++ b/src/ifc2ca/templates/codeaster/finish_process.py @@ -0,0 +1,4 @@ +# STEP: CONCLUDE STUDY +# code_aster.close() +FIN() +{{"\n"}} diff --git a/src/ifc2ca/templates/codeaster/read_mesh.py b/src/ifc2ca/templates/codeaster/read_mesh.py new file mode 100644 index 0000000000..0799c24c89 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/read_mesh.py @@ -0,0 +1,3 @@ +# STEP: READ MED FILE +mesh = LIRE_MAILLAGE(FORMAT="MED", UNITE=20) +{{"\n"}} diff --git a/src/ifc2ca/templates/codeaster/run_analysis.py b/src/ifc2ca/templates/codeaster/run_analysis.py new file mode 100644 index 0000000000..4f1ca8a066 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/run_analysis.py @@ -0,0 +1,27 @@ +# STEP: RUN ANALYSIS +{{ res_Bld }} = MECA_STATIQUE( + MODELE = model, + CHAM_MATER = material, + CARA_ELEM = element, + LIST_INST = {{ analysis_time }}, + EXCIT = ( + _F( + CHARGE = connection + ), + _F( + CHARGE = {{ load }} + ) + ), + SOLVEUR=_F( + NPREC=12, + RESI_RELA=1e-1, + STOP_SINGULIER='NON', + ) +) + +{{ res_Bld }} = CALC_CHAMP( + reuse = {{ res_Bld }}, + RESULTAT = {{ res_Bld }}, + CONTRAINTE=('EFGE_NOEU', ), +) +{{ "\n" }} diff --git a/src/ifc2ca/templates/codeaster/start_process.py b/src/ifc2ca/templates/codeaster/start_process.py new file mode 100644 index 0000000000..dc7a883a54 --- /dev/null +++ b/src/ifc2ca/templates/codeaster/start_process.py @@ -0,0 +1,24 @@ +# Ifc2CA - IFC Code_Aster utility +# Copyright (C) 2020, 2021, 2023, 2024 Ioannis P. Christovasilis +# +# This file is part of Ifc2CA. +# +# Ifc2CA is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Ifc2CA 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Ifc2CA. If not, see . + + +# STEP: INITIALIZE STUDY +# import code_aster +# code_aster.init() +DEBUT() +{{"\n"}} diff --git a/src/ifc2ca/scriptSalome.py b/src/ifc2ca/templates/salome/scriptSalome.py similarity index 57% rename from src/ifc2ca/scriptSalome.py rename to src/ifc2ca/templates/salome/scriptSalome.py index 35e295760f..756bbcebdd 100644 --- a/src/ifc2ca/scriptSalome.py +++ b/src/ifc2ca/templates/salome/scriptSalome.py @@ -1,5 +1,5 @@ # Ifc2CA - IFC Code_Aster utility -# Copyright (C) 2020, 2021 Ioannis P. Christovasilis +# Copyright (C) 2020, 2021, 2023, 2024 Ioannis P. Christovasilis # # This file is part of Ifc2CA. # @@ -16,26 +16,31 @@ # You should have received a copy of the GNU Lesser General Public License # along with Ifc2CA. If not, see . -from __future__ import division -from __future__ import print_function +import itertools +import json import os import time -import json +from pathlib import Path + +import numpy as np import salome import salome_notebook import salome_version -import numpy as np -import itertools -from pathlib import Path flatten = itertools.chain.from_iterable +mesh_size = {{ mesh_size }} +med_path = r"{{ med_path }}" +json_path = r"{{ json_path }}" + +with open(json_path, "r") as f: + data = json.load(f) + class MODEL: - def __init__(self, dataFilename, medFilename, meshSize): - self.dataFilename = dataFilename - self.medFilename = medFilename - self.meshSize = meshSize + def __init__(self): + self.medFilename = med_path + self.mesh_size = mesh_size self.tolLoc = 0 self.mesh = None self.meshNodes = None @@ -82,29 +87,22 @@ class MODEL: return self.geompy.MakeFaceWires(LineList, 1) - def makeObject(self, geometry, geometryType): - if geometryType == "point": + def makeObject(self, geometry, geometry_type): + if geometry_type == "Vertex": return self.makePoint(geometry) - if geometryType == "line": + if geometry_type == "Edge": return self.makeLine(geometry) - if geometryType == "surface": + if geometry_type == "Face": return self.makeFace(geometry) - def makePartition(self, objects, geometryType): - if geometryType == "point": + def makePartition(self, objects, geometry_type): + if geometry_type == "Vertex": shapeType = "VERTEX" - if geometryType == "line": + if geometry_type == "Edge": shapeType = "EDGE" - if geometryType == "surface": + if geometry_type == "Face": shapeType = "FACE" - return self.geompy.MakePartition( - objects, [], [], [], self.geompy.ShapeType[shapeType], 0, [], 1 - ) - - def getLinkGeometry(self, ecc, orientation, finalPoint): - vector = np.array(orientation).transpose().dot(ecc["vector"]) - initialPoint = (np.array(finalPoint) - vector).tolist() - return [initialPoint, finalPoint] + return self.geompy.MakePartition(objects, [], [], [], self.geompy.ShapeType[shapeType], 0, [], 1) def length(self, geometry): return ( @@ -115,18 +113,17 @@ class MODEL: def create(self): # Read data from input file - with open(self.dataFilename) as dataFile: - data = json.load(dataFile) + # data = data - elements = data["elements"] - connections = data["connections"] + self.elements = elements = data["elements"] + self.connections = connections = data["connections"] # --> Delete this reference data and repopulate it with the objects # while going through elements for conn in connections: - conn["relatedElements"] = [] + conn["related_elements"] = [] # End <-- - meshSize = self.meshSize + mesh_size = self.mesh_size dec = 7 # 4 decimals for length in mm tol = 10 ** (-dec - 3 + 1) @@ -134,7 +131,7 @@ class MODEL: self.tolLoc = tol * 10 * 2 tolLoc = self.tolLoc - NEW_SALOME = int(salome_version.getVersion()[0]) >= 9 + self.NEW_SALOME = NEW_SALOME = int(salome_version.getVersion()[0]) >= 9 salome.salome_init() theStudy = salome.myStudy notebook = salome_notebook.NoteBook(theStudy) @@ -142,10 +139,11 @@ class MODEL: ### ### GEOM component ### - import GEOM - from salome.geom import geomBuilder import math + + import GEOM import SALOMEDS + from salome.geom import geomBuilder gg = salome.ImportComponentGUI("GEOM") if NEW_SALOME: @@ -163,9 +161,9 @@ class MODEL: geompy.addToStudy(OY, "OY") geompy.addToStudy(OZ, "OZ") - if len([e for e in elements if e["geometryType"] == "line"]) > 0: + if len([e for e in elements if e["geometry_type"] == "Edge"]) > 0: buildingShapeType = "EDGE" - if len([e for e in elements if e["geometryType"] == "surface"]) > 0: + if len([e for e in elements if e["geometry_type"] == "Face"]) > 0: buildingShapeType = "FACE" ### Define entities ### @@ -175,31 +173,23 @@ class MODEL: # Loop 1 for el in elements: - el["elemObj"] = self.makeObject(el["geometry"], el["geometryType"]) + el["elemObj"] = self.makeObject(el["geometry"], el["geometry_type"]) el["connObjs"] = [None for _ in el["connections"]] el["linkObjs"] = [None for _ in el["connections"]] el["linkPointObjs"] = [[None, None] for _ in el["connections"]] for j, rel in enumerate(el["connections"]): - conn = [ - c for c in connections if c["referenceName"] == rel["relatedConnection"] - ][0] + conn = [c for c in connections if c["ref_id"] == rel["related_connection"]][0] if rel["eccentricity"]: - rel["index"] = len(conn["relatedElements"]) + 1 - conn["relatedElements"].append(rel) + rel["index"] = len(conn["related_elements"]) + 1 + conn["related_elements"].append(rel) if not rel["eccentricity"]: - el["connObjs"][j] = self.makeObject( - conn["geometry"], conn["geometryType"] - ) + el["connObjs"][j] = self.makeObject(conn["geometry"], conn["geometry_type"]) else: - if conn["geometryType"] == "point": - geometry = self.getLinkGeometry( - rel["eccentricity"], el["orientation"], conn["geometry"] - ) - el["connObjs"][j] = self.makeObject( - geometry[0], conn["geometryType"] - ) + if conn["geometry_type"] == "Vertex": + geometry = rel["eccentricity"]["point_on_element"], conn["geometry"] + el["connObjs"][j] = self.makeObject(geometry[0], conn["geometry_type"]) el["linkPointObjs"][j][0] = self.geompy.MakeVertex( geometry[0][0], geometry[0][1], geometry[0][2] @@ -211,27 +201,18 @@ class MODEL: el["linkPointObjs"][j][0], el["linkPointObjs"][j][1] ) else: - print( - "Eccentricity defined for a %s geometryType" - % conn["geometryType"] - ) - el["partObj"] = self.makePartition( - [el["elemObj"]] + el["connObjs"], el["geometryType"] - ) + print("Eccentricity defined for a %s geometry_type" % conn["geometry_type"]) + el["partObj"] = self.makePartition([el["elemObj"]] + el["connObjs"], el["geometry_type"]) el["elemObj"] = geompy.GetInPlace(el["partObj"], el["elemObj"], True) for j, rel in enumerate(el["connections"]): - el["connObjs"][j] = geompy.GetInPlace( - el["partObj"], el["connObjs"][j], True - ) + el["connObjs"][j] = geompy.GetInPlace(el["partObj"], el["connObjs"][j], True) for conn in connections: - conn["connObj"] = self.makeObject(conn["geometry"], conn["geometryType"]) + conn["connObj"] = self.makeObject(conn["geometry"], conn["geometry_type"]) # Make assemble of Building Object bldObjs = [] bldObjs.extend([el["partObj"] for el in elements]) - bldObjs.extend( - flatten([[link for link in el["linkObjs"] if link] for el in elements]) - ) + bldObjs.extend(flatten([[link for link in el["linkObjs"] if link] for el in elements])) bldObjs.extend([conn["connObj"] for conn in connections]) bldComp = geompy.MakeCompound(bldObjs) @@ -240,59 +221,55 @@ class MODEL: # Loop 2 for el in elements: - # geompy.addToStudy(el['partObj'], self.getGroupName(el['referenceName'])) - geompy.addToStudyInFather( - el["partObj"], el["elemObj"], self.getGroupName(el["referenceName"]) - ) + # geompy.addToStudy(el['partObj'], self.getGroupName(el['ref_id'])) + geompy.addToStudyInFather(el["partObj"], el["elemObj"], self.getGroupName(el["ref_id"])) for j, rel in enumerate(el["connections"]): - conn = [ - c for c in connections if c["referenceName"] == rel["relatedConnection"] - ][0] + conn = [c for c in connections if c["ref_id"] == rel["related_connection"]][0] rel["conn_string"] = None - if conn["geometryType"] == "point": + if conn["geometry_type"] == "Vertex": rel["conn_string"] = "_0DC_" - if conn["geometryType"] == "line": + if conn["geometry_type"] == "Edge": rel["conn_string"] = "_1DC_" - if conn["geometryType"] == "surface": + if conn["geometry_type"] == "Face": rel["conn_string"] = "_2DC_" geompy.addToStudyInFather( el["partObj"], el["connObjs"][j], - self.getGroupName(el["referenceName"]) - + rel["conn_string"] - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + rel["conn_string"] + self.getGroupName(rel["related_connection"]), ) if rel["eccentricity"]: pass - # geompy.addToStudy(el['linkObjs'][j], self.getGroupName(el['referenceName']) + '_1DR_' + self.getGroupName(rel['relatedConnection'])) - # geompy.addToStudyInFather(el['linkObjs'][j], el['linkPointObjs'][j][0], self.getGroupName(rel['relatedConnection']) + '_0DC_' + self.getGroupName(el['referenceName'])) - # geompy.addToStudyInFather(el['linkObjs'][j], el['linkPointObjs'][j][0], self.getGroupName(rel['relatedConnection']) + '_0DC_%g' % rel['index']) + # geompy.addToStudy(el['linkObjs'][j], self.getGroupName(el['ref_id']) + '_1DR_' + self.getGroupName(rel['related_connection'])) + # geompy.addToStudyInFather(el['linkObjs'][j], el['linkPointObjs'][j][0], self.getGroupName(rel['related_connection']) + '_0DC_' + self.getGroupName(el['ref_id'])) + # geompy.addToStudyInFather(el['linkObjs'][j], el['linkPointObjs'][j][0], self.getGroupName(rel['related_connection']) + '_0DC_%g' % rel['index']) for conn in connections: - # geompy.addToStudy(conn['connObj'], self.getGroupName(conn['referenceName'])) - geompy.addToStudyInFather( - conn["connObj"], conn["connObj"], self.getGroupName(conn["referenceName"]) - ) + # geompy.addToStudy(conn['connObj'], self.getGroupName(conn['ref_id'])) + geompy.addToStudyInFather(conn["connObj"], conn["connObj"], self.getGroupName(conn["ref_id"])) elapsed_time = time.time() - init_time init_time += elapsed_time print("Building Geometry Defined in %g sec" % (elapsed_time)) # Define and add groups for all curve and surface members - if len([e for e in elements if e["geometryType"] == "line"]) > 0: + if len([e for e in elements if e["geometry_type"] == "Edge"]) > 0: # Make compound of requested group - compoundTemp = geompy.MakeCompound( - [e["elemObj"] for e in elements if e["geometryType"] == "line"] - ) + compoundTemp = geompy.MakeCompound([e["elemObj"] for e in elements if e["geometry_type"] == "Edge"]) # Define group object and add to study curveCompound = geompy.GetInPlace(bldComp, compoundTemp, True) geompy.addToStudyInFather(bldComp, curveCompound, "CurveMembers") - if len([e for e in elements if e["geometryType"] == "surface"]) > 0: + rigid_links = list(flatten([[link for link in el["linkObjs"] if link] for el in elements])) + if len(rigid_links) > 0: # Make compound of requested group - compoundTemp = geompy.MakeCompound( - [e["elemObj"] for e in elements if e["geometryType"] == "surface"] - ) + compoundTemp = geompy.MakeCompound(rigid_links) + # Define group object and add to study + rigidLinkCompound = geompy.GetInPlace(bldComp, compoundTemp, True) + geompy.addToStudyInFather(bldComp, rigidLinkCompound, "RigidLinks") + + if len([e for e in elements if e["geometry_type"] == "Face"]) > 0: + # Make compound of requested group + compoundTemp = geompy.MakeCompound([e["elemObj"] for e in elements if e["geometry_type"] == "Face"]) # Define group object and add to study surfaceCompound = geompy.GetInPlace(bldComp, compoundTemp, True) geompy.addToStudyInFather(bldComp, surfaceCompound, "SurfaceMembers") @@ -300,45 +277,34 @@ class MODEL: # Loop 3 for el in elements: # el['partObj'] = geompy.RestoreGivenSubShapes(bldComp, [el['partObj']], GEOM.FSM_GetInPlace, False, False)[0] - geompy.addToStudyInFather( - bldComp, el["elemObj"], self.getGroupName(el["referenceName"]) - ) + geompy.addToStudyInFather(bldComp, el["elemObj"], self.getGroupName(el["ref_id"])) for j, rel in enumerate(el["connections"]): geompy.addToStudyInFather( bldComp, el["connObjs"][j], - self.getGroupName(el["referenceName"]) - + rel["conn_string"] - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + rel["conn_string"] + self.getGroupName(rel["related_connection"]), ) if rel["eccentricity"]: # point geometry geompy.addToStudyInFather( bldComp, el["linkObjs"][j], - self.getGroupName(el["referenceName"]) - + "_1DR_" - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + "_1DR_" + self.getGroupName(rel["related_connection"]), ) geompy.addToStudyInFather( bldComp, el["linkPointObjs"][j][0], - self.getGroupName(rel["relatedConnection"]) - + "_0DC_" - + self.getGroupName(el["referenceName"]), + self.getGroupName(rel["related_connection"]) + "_0DC_" + self.getGroupName(el["ref_id"]), ) geompy.addToStudyInFather( bldComp, el["linkPointObjs"][j][1], - self.getGroupName(rel["relatedConnection"]) - + "_0DC_%g" % rel["index"], + self.getGroupName(rel["related_connection"]) + "_0DC_%g" % rel["index"], ) for conn in connections: # conn['connObj'] = geompy.RestoreGivenSubShapes(bldComp, [conn['connObj']], GEOM.FSM_GetInPlace, False, False)[0] - geompy.addToStudyInFather( - bldComp, conn["connObj"], self.getGroupName(conn["referenceName"]) - ) + geompy.addToStudyInFather(bldComp, conn["connObj"], self.getGroupName(conn["ref_id"])) elapsed_time = time.time() - init_time init_time += elapsed_time @@ -359,15 +325,15 @@ class MODEL: smesh = smeshBuilder.New(theStudy) bldMesh = smesh.Mesh(bldComp) Regular_1D = bldMesh.Segment() - Local_Length_1 = Regular_1D.LocalLength(meshSize, None, tolLoc) + Local_Length_1 = Regular_1D.LocalLength(mesh_size, None, tolLoc) if buildingShapeType == "FACE": NETGEN2D_ONLY = bldMesh.Triangle(algo=smeshBuilder.NETGEN_2D) NETGEN2D_Pars = NETGEN2D_ONLY.Parameters() - NETGEN2D_Pars.SetMaxSize(meshSize) + NETGEN2D_Pars.SetMaxSize(mesh_size) NETGEN2D_Pars.SetOptimize(1) NETGEN2D_Pars.SetFineness(2) - NETGEN2D_Pars.SetMinSize(meshSize / 5.0) + NETGEN2D_Pars.SetMinSize(mesh_size / 5.0) NETGEN2D_Pars.SetUseSurfaceCurvature(1) NETGEN2D_Pars.SetQuadAllowed(1) NETGEN2D_Pars.SetSecondOrder(0) @@ -383,142 +349,111 @@ class MODEL: smesh.SetName(NETGEN2D_ONLY.GetAlgorithm(), "NETGEN2D_ONLY") smesh.SetName(NETGEN2D_Pars, "NETGEN2D_Pars") - smesh.SetName(bldMesh.GetMesh(), "bldMesh") + smesh.SetName(bldMesh.GetMesh(), "{{ mesh_name }}") elapsed_time = time.time() - init_time init_time += elapsed_time print("Meshing Operations Completed in %g sec" % (elapsed_time)) # Define and add groups for all curve and surface members - if len([e for e in elements if e["geometryType"] == "line"]) > 0: + if len([e for e in elements if e["geometry_type"] == "Edge"]) > 0: tempgroup = bldMesh.GroupOnGeom(curveCompound, "CurveMembers", SMESH.EDGE) smesh.SetName(tempgroup, "CurveMembers") - if len([e for e in elements if e["geometryType"] == "surface"]) > 0: + if len(rigid_links) > 0: + tempgroup = bldMesh.GroupOnGeom(rigidLinkCompound, "RigidLinks", SMESH.EDGE) + smesh.SetName(tempgroup, "RigidLinks") + + if len([e for e in elements if e["geometry_type"] == "Face"]) > 0: tempgroup = bldMesh.GroupOnGeom(surfaceCompound, "SurfaceMembers", SMESH.FACE) smesh.SetName(tempgroup, "SurfaceMembers") # Define groups in Mesh for el in elements: - if el["geometryType"] == "line": + if el["geometry_type"] == "Edge": shapeType = SMESH.EDGE - if el["geometryType"] == "surface": + if el["geometry_type"] == "Face": shapeType = SMESH.FACE - tempgroup = bldMesh.GroupOnGeom( - el["elemObj"], self.getGroupName(el["referenceName"]), shapeType - ) - smesh.SetName(tempgroup, self.getGroupName(el["referenceName"])) + tempgroup = bldMesh.GroupOnGeom(el["elemObj"], self.getGroupName(el["ref_id"]), shapeType) + smesh.SetName(tempgroup, self.getGroupName(el["ref_id"])) + # tempgroup = bldMesh.GroupOnGeom(el["elemObj"], self.getGroupName(el["ref_id"]), SMESH.NODE) + # smesh.SetName(tempgroup, self.getGroupName(el["ref_id"])) for j, rel in enumerate(el["connections"]): tempgroup = bldMesh.GroupOnGeom( el["connObjs"][j], - self.getGroupName(el["referenceName"]) - + rel["conn_string"] - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + rel["conn_string"] + self.getGroupName(rel["related_connection"]), SMESH.NODE, ) smesh.SetName( tempgroup, - self.getGroupName(el["referenceName"]) - + rel["conn_string"] - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + rel["conn_string"] + self.getGroupName(rel["related_connection"]), ) - rel["node"] = ( - bldMesh.GetIDSource(tempgroup.GetNodeIDs(), SMESH.NODE) - ).GetIDs()[0] + rel["node"] = (bldMesh.GetIDSource(tempgroup.GetNodeIDs(), SMESH.NODE)).GetIDs()[0] if rel["eccentricity"]: tempgroup = bldMesh.GroupOnGeom( el["linkObjs"][j], - self.getGroupName(el["referenceName"]) - + "_1DR_" - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + "_1DR_" + self.getGroupName(rel["related_connection"]), SMESH.EDGE, ) smesh.SetName( tempgroup, - self.getGroupName(el["referenceName"]) - + "_1DR_" - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + "_1DR_" + self.getGroupName(rel["related_connection"]), ) tempgroup = bldMesh.GroupOnGeom( el["linkPointObjs"][j][0], - self.getGroupName(rel["relatedConnection"]) - + "_0DC_" - + self.getGroupName(el["referenceName"]), + self.getGroupName(rel["related_connection"]) + "_0DC_" + self.getGroupName(el["ref_id"]), SMESH.NODE, ) smesh.SetName( tempgroup, - self.getGroupName(rel["relatedConnection"]) - + "_0DC_" - + self.getGroupName(el["referenceName"]), + self.getGroupName(rel["related_connection"]) + "_0DC_" + self.getGroupName(el["ref_id"]), ) - rel["eccNode"] = ( - bldMesh.GetIDSource(tempgroup.GetNodeIDs(), SMESH.NODE) - ).GetIDs()[0] + rel["eccNode"] = (bldMesh.GetIDSource(tempgroup.GetNodeIDs(), SMESH.NODE)).GetIDs()[0] tempgroup = bldMesh.GroupOnGeom( el["linkPointObjs"][j][1], - self.getGroupName(rel["relatedConnection"]) + self.getGroupName(rel["related_connection"]) + "_0DC_" - + self.getGroupName(rel["relatedConnection"]), + + self.getGroupName(rel["related_connection"]), SMESH.NODE, ) smesh.SetName( tempgroup, - self.getGroupName(rel["relatedConnection"]) - + "_0DC_%g" % rel["index"], + self.getGroupName(rel["related_connection"]) + "_0DC_%g" % rel["index"], ) for conn in connections: - tempgroup = bldMesh.GroupOnGeom( - conn["connObj"], self.getGroupName(conn["referenceName"]), SMESH.NODE - ) - smesh.SetName(tempgroup, self.getGroupName(conn["referenceName"])) + tempgroup = bldMesh.GroupOnGeom(conn["connObj"], self.getGroupName(conn["ref_id"]), SMESH.NODE) + smesh.SetName(tempgroup, self.getGroupName(conn["ref_id"])) nodesId = bldMesh.GetIDSource(tempgroup.GetNodeIDs(), SMESH.NODE) - tempgroup = bldMesh.Add0DElementsToAllNodes( - nodesId, self.getGroupName(conn["referenceName"]) - ) - smesh.SetName(tempgroup, self.getGroupName(conn["referenceName"] + "_0D")) - if conn["geometryType"] == "point": + tempgroup = bldMesh.Add0DElementsToAllNodes(nodesId, self.getGroupName(conn["ref_id"])) + smesh.SetName(tempgroup, self.getGroupName(conn["ref_id"] + "_0D")) + if conn["geometry_type"] == "Vertex": conn["node"] = nodesId.GetIDs()[0] - if conn["geometryType"] == "line": - tempgroup = bldMesh.GroupOnGeom( - conn["connObj"], self.getGroupName(conn["referenceName"]), SMESH.EDGE - ) - smesh.SetName(tempgroup, self.getGroupName(conn["referenceName"])) - if conn["geometryType"] == "surface": - tempgroup = bldMesh.GroupOnGeom( - conn["connObj"], self.getGroupName(conn["referenceName"]), SMESH.FACE - ) - smesh.SetName(tempgroup, self.getGroupName(conn["referenceName"])) + if conn["geometry_type"] == "Edge": + tempgroup = bldMesh.GroupOnGeom(conn["connObj"], self.getGroupName(conn["ref_id"]), SMESH.EDGE) + smesh.SetName(tempgroup, self.getGroupName(conn["ref_id"])) + if conn["geometry_type"] == "Face": + tempgroup = bldMesh.GroupOnGeom(conn["connObj"], self.getGroupName(conn["ref_id"]), SMESH.FACE) + smesh.SetName(tempgroup, self.getGroupName(conn["ref_id"])) # create 1D SEG2 spring elements for el in elements: for j, rel in enumerate(el["connections"]): - conn = [ - c for c in connections if c["referenceName"] == rel["relatedConnection"] - ][0] - if conn["geometryType"] == "point": + conn = [c for c in connections if c["ref_id"] == rel["related_connection"]][0] + if conn["geometry_type"] == "Vertex": grpName = bldMesh.CreateEmptyGroup( SMESH.EDGE, - self.getGroupName(el["referenceName"]) - + "_1DS_" - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + "_1DS_" + self.getGroupName(rel["related_connection"]), ) smesh.SetName( grpName, - self.getGroupName(el["referenceName"]) - + "_1DS_" - + self.getGroupName(rel["relatedConnection"]), + self.getGroupName(el["ref_id"]) + "_1DS_" + self.getGroupName(rel["related_connection"]), ) if not rel["eccentricity"]: - conn = [ - conn - for conn in connections - if conn["referenceName"] == rel["relatedConnection"] - ][0] + conn = [conn for conn in connections if conn["ref_id"] == rel["related_connection"]][0] grpName.Add([bldMesh.AddEdge([conn["node"], rel["node"]])]) else: grpName.Add([bldMesh.AddEdge([rel["eccNode"], rel["node"]])]) @@ -545,26 +480,30 @@ class MODEL: except: print("ExportMED() failed. Invalid file name?") - if salome.sg.hasDesktop(): - if NEW_SALOME: - salome.sg.updateObjBrowser() - else: - salome.sg.updateObjBrowser(1) + # if salome.sg.hasDesktop(): + # if NEW_SALOME: + # salome.sg.updateObjBrowser() + # else: + # salome.sg.updateObjBrowser(1) elapsed_time = init_time - start_time print("ALL Operations Completed in %g sec" % (elapsed_time)) -if __name__ == "__main__": - fileNames = ["structure_01"] - files = fileNames +model = MODEL() - meshSize = 0.1 +for el in model.elements: + for j, conn in enumerate(el["connections"]): + d = model.geompy.MinDistance(el["elemObj"], el["connObjs"][j]) + if d > 0: + print(f'NOTE: Element {el["ref_id"]} and connection {conn["ref_id"]} have a distance of {d}') + # elif d == 0: + # print( + # f'SUCCESS: Element {el["ref_id"]} and connection {conn["ref_id"]} have a distance of {d}' + # ) - for fileName in files: - BASE_PATH = Path( - "/home/jesusbill/Dev-Projects/github.com/IfcOpenShell/analysis-models/models/" - ) - DATAFILENAME = BASE_PATH / fileName / f"{fileName}.json" - MEDFILENAME = BASE_PATH / fileName / f"{fileName}.med" - model = MODEL(DATAFILENAME, str(MEDFILENAME), meshSize) +if salome.sg.hasDesktop(): + if model.NEW_SALOME: + salome.sg.updateObjBrowser() + else: + salome.sg.updateObjBrowser(1)