add "bonded" case for mesh and run

src/ifc2ca/ifc2ca.py

@@ -17,7 +17,9 @@ class IFC2CA:
         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")
+            connections = self.get_structural_items(
+                model, item_type="IfcStructuralConnection"
+            )
 
             materialdb = []
             materials = list(dict.fromkeys([e["material"] for e in elements]))
@@ -25,16 +27,24 @@ class IFC2CA:
                 id = int(mat.split("|")[1])
                 material = self.get_material_properties(self.file.by_id(id))
                 material["relatedElements"] = [
-                    e["ifcName"] for e in elements if "material" in e and e["material"] == mat
+                    e["ifcName"]
+                    for e in elements
+                    if "material" in e and e["material"] == mat
                 ]
                 materialdb.append(material)
 
             profiledb = []
-            profiles = list(dict.fromkeys([e["profile"] for e in elements if "profile" in e]))
+            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["ifcName"] for e in elements if "profile" in e and e["profile"] == prof]
+                profile["relatedElements"] = [
+                    e["ifcName"]
+                    for e in elements
+                    if "profile" in e and e["profile"] == prof
+                ]
                 profiledb.append(profile)
 
             self.result = {
@@ -47,11 +57,11 @@ class IFC2CA:
                 "warnings": self.warnings,
             }
 
-            print('Model "%s" converted' % model.Name)
+            print(f"Model {model.Name} converted")
-            print("Number of elements: ", len(elements))
+            print(f"Number of elements: {len(elements)}")
-            print("Number of connections: ", len(connections))
+            print(f"Number of connections: {len(connections)}")
-            print("Number of materials: ", len(materialdb))
+            print(f"Number of materials: {len(materialdb)}")
-            print("Number of profiles: ", len(profiledb))
+            print(f"Number of profiles: {len(profiledb)}")
             print("")
 
             break
@@ -75,12 +85,19 @@ class IFC2CA:
             material_profile = self.get_material_profile(item)
             if not representation:
                 self.warnings.append(
-                    "No representation defined for %s. Member excluded" % (item.is_a() + "|" + str(item.id()))
+                    "No representation defined for %s. Member excluded"
+                    % (item.is_a() + "|" + str(item.id()))
                 )
                 return
             if not material_profile:
-                self.warnings.append("No material defined for in %s" % (item.is_a() + "|" + str(item.id())))
+                self.warnings.append(
-                self.warnings.append("No profile defined for in %s" % (item.is_a() + "|" + str(item.id())))
+                    "No material defined for in %s"
+                    % (item.is_a() + "|" + str(item.id()))
+                )
+                self.warnings.append(
+                    "No profile defined for in %s"
+                    % (item.is_a() + "|" + str(item.id()))
+                )
                 materialId = None
                 profileId = None
             else:
@@ -99,21 +116,32 @@ class IFC2CA:
             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:
+                    if (
-                        print(np.linalg.norm(np.array(c["eccentricity"]["pointOnElement"])), ">", length)
+                        np.linalg.norm(np.array(c["eccentricity"]["pointOnElement"]))
-                        self.warnings.append("Eccentricity in %s corrected" % (item.is_a() + "|" + str(item.id())))
+                        > length + self.tol
+                    ):
+                        print(
+                            f"{np.linalg.norm(np.array(c['eccentricity']['pointOnElement']))} > {length}"
+                        )
+                        self.warnings.append(
+                            f"Eccentricity in {item.is_a()}|{str(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)
+                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
+                            c["eccentricity"]["vector"],
+                            transformation,
+                            include_translation=False,
                         )
 
             return {
@@ -134,11 +162,15 @@ class IFC2CA:
             material = self.get_material_profile(item)
             if not representation:
                 self.warnings.append(
-                    "No representation defined for %s. Member excluded" % (item.is_a() + "|" + str(item.id()))
+                    "No representation defined for %s. Member excluded"
+                    % (item.is_a() + "|" + str(item.id()))
                 )
                 return
             if not material:
-                self.warnings.append("No material defined for in %s" % (item.is_a() + "|" + str(item.id())))
+                self.warnings.append(
+                    "No material defined for in %s"
+                    % (item.is_a() + "|" + str(item.id()))
+                )
                 materialId = None
             else:
                 materialId = material.is_a() + "|" + str(material.id())
@@ -151,7 +183,9 @@ class IFC2CA:
                     conn["orientation"] = orientation
             if transformation:
                 geometry = self.transform_vectors(geometry, transformation)
-                orientation = self.transform_vectors(orientation, transformation, include_translation=False)
+                orientation = self.transform_vectors(
+                    orientation, transformation, include_translation=False
+                )
                 for c in connections:
                     c["orientation"] = self.transform_vectors(
                         c["orientation"], transformation, include_translation=False
@@ -174,7 +208,8 @@ class IFC2CA:
             representation = self.get_representation(item, "Vertex")
             if not representation:
                 self.warnings.append(
-                    "No representation defined for %s. Connection excluded" % (item.is_a() + "|" + str(item.id()))
+                    "No representation defined for %s. Connection excluded"
+                    % (item.is_a() + "|" + str(item.id()))
                 )
                 return
 
@@ -184,7 +219,9 @@ class IFC2CA:
                 orientation = np.eye(3).tolist()
             if transformation:
                 geometry = self.transform_vectors(geometry, transformation)
-                orientation = self.transform_vectors(orientation, transformation, include_translation=False)
+                orientation = self.transform_vectors(
+                    orientation, transformation, include_translation=False
+                )
 
             return {
                 "ifcName": item.is_a() + "|" + str(item.id()),
@@ -194,14 +231,18 @@ class IFC2CA:
                 "geometry": geometry,
                 "orientation": orientation,
                 "appliedCondition": self.get_connection_input(item, "point"),
-                "relatedElements": [con.is_a() + "|" + str(con.id()) for con in item.ConnectsStructuralMembers],
+                "relatedElements": [
+                    con.is_a() + "|" + str(con.id())
+                    for con in item.ConnectsStructuralMembers
+                ],
             }
 
         elif item.is_a("IfcStructuralCurveConnection"):
             representation = self.get_representation(item, "Edge")
             if not representation:
                 self.warnings.append(
-                    "No representation defined for %s. Connection excluded" % (item.is_a() + "|" + str(item.id()))
+                    "No representation defined for %s. Connection excluded"
+                    % (item.is_a() + "|" + str(item.id()))
                 )
                 return
 
@@ -211,7 +252,9 @@ class IFC2CA:
                 orientation = np.eye(3).tolist()
             if transformation:
                 geometry = self.transform_vectors(geometry, transformation)
-                orientation = self.transform_vectors(orientation, transformation, include_translation=False)
+                orientation = self.transform_vectors(
+                    orientation, transformation, include_translation=False
+                )
 
             return {
                 "ifcName": item.is_a() + "|" + str(item.id()),
@@ -221,7 +264,10 @@ class IFC2CA:
                 "geometry": geometry,
                 "orientation": orientation,
                 "appliedCondition": self.get_connection_input(item, "line"),
-                "relatedElements": [con.is_a() + "|" + str(con.id()) for con in item.ConnectsStructuralMembers],
+                "relatedElements": [
+                    con.is_a() + "|" + str(con.id())
+                    for con in item.ConnectsStructuralMembers
+                ],
             }
 
     def get_transformation(self, placement):
@@ -229,11 +275,15 @@ class IFC2CA:
             return None
         if placement.is_a("IfcLocalPlacement"):
             if placement.PlacementRelTo:
-                print("Warning! Object Placement with PlacementRelTo attribute is not supported and will be neglected")
+                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)
+                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)
@@ -253,10 +303,13 @@ class IFC2CA:
                 and np.allclose(zAxis, np.array([0.0, 0.0, 1.0]))
             ):
                 return None
-            return {"location": location, "rotationMatrix": np.array([xAxis, yAxis, zAxis]).transpose()}
+            return {
+                "location": location,
+                "rotationMatrix": np.array([xAxis, yAxis, zAxis]).transpose(),
+            }
         else:
             print(
-                "Warning! Object Placement is of type %s, which is not supported. Default considered" % placement.is_a()
+                f"Warning! Object Placement is of type {placement.is_a()}, which is not supported. Default considered"
             )
             return None
 
@@ -264,11 +317,13 @@ class IFC2CA:
         if not element.Representation:
             return None
         for representation in element.Representation.Representations:
-            rep = self.get_specific_representation(representation, "Reference", rep_type)
+            rep = self.get_specific_representation(
+                representation, "Reference", rep_type
+            )
             if rep:
                 return rep
         else:
-            # print('Trying without rep identifier')
+            # print("Trying without rep identifier")
             for representation in element.Representation.Representations:
                 rep = self.get_specific_representation(representation, None, rep_type)
                 if rep:
@@ -281,7 +336,9 @@ class IFC2CA:
             return representation
         if representation.RepresentationType == "MappedRepresentation":
             return self.get_specific_representation(
-                representation.Items[0].MappingSource.MappedRepresentation, rep_id, rep_type
+                representation.Items[0].MappingSource.MappedRepresentation,
+                rep_id,
+                rep_type,
             )
 
     def get_geometry(self, representation):
@@ -299,7 +356,9 @@ class IFC2CA:
             edges = item.Bounds[0].Bound.EdgeList
             coords = []
             for edge in edges:
-                coords.append(self.get_coordinate(edge.EdgeElement.EdgeStart.VertexGeometry))
+                coords.append(
+                    self.get_coordinate(edge.EdgeElement.EdgeStart.VertexGeometry)
+                )
             return coords
 
         elif item.is_a("IfcVertexPoint"):
@@ -328,7 +387,9 @@ class IFC2CA:
     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)
+        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)
@@ -347,7 +408,9 @@ class IFC2CA:
             return orientation
 
     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
+        if not any(
+            isinstance(el, list) for el in geometry
+        ):  # single point which contains no list
             geometry = [geometry]
         globalGeometry = []
 
@@ -453,13 +516,18 @@ class IFC2CA:
             {
                 "ifcName": rel.is_a() + "|" + str(rel.id()),
                 "id": rel.GlobalId,
-                "relatingElement": rel.RelatingStructuralMember.is_a() + "|" + str(rel.RelatingStructuralMember.id()),
+                "relatingElement": rel.RelatingStructuralMember.is_a()
+                + "|"
+                + str(rel.RelatingStructuralMember.id()),
                 "relatedConnection": rel.RelatedStructuralConnection.is_a()
                 + "|"
                 + str(rel.RelatedStructuralConnection.id()),
                 "orientation": self.get_0D_orientation(rel.ConditionCoordinateSystem),
                 "appliedCondition": self.get_connection_input(
-                    rel, self.get_geometry_type_from_connection(rel.RelatedStructuralConnection)
+                    rel,
+                    self.get_geometry_type_from_connection(
+                        rel.RelatedStructuralConnection
+                    ),
                 ),
                 "eccentricity": None
                 if not rel.is_a("IfcRelConnectsWithEccentricity")
@@ -475,7 +543,9 @@ class IFC2CA:
                         if not rel.ConnectionConstraint.EccentricityInZ
                         else rel.ConnectionConstraint.EccentricityInZ,
                     ],
-                    "pointOnElement": self.get_coordinate(rel.ConnectionConstraint.PointOnRelatingElement),
+                    "pointOnElement": self.get_coordinate(
+                        rel.ConnectionConstraint.PointOnRelatingElement
+                    ),
                 },
             }
             for rel in itemList
@@ -532,11 +602,23 @@ class IFC2CA:
             Iz = (2 * tf) * (b ** 3) / 12 + (h - 2 * tf) * (tw ** 3) / 12
             Jx = 1 / 3 * ((h - tf) * (tw ** 3) + 2 * b * (tf ** 3))
 
-            return {"crossSectionArea": A, "momentOfInertiaY": Iy, "momentOfInertiaZ": Iz, "torsionalConstantX": Jx}
+            return {
+                "crossSectionArea": A,
+                "momentOfInertiaY": Iy,
+                "momentOfInertiaZ": Iz,
+                "torsionalConstantX": Jx,
+            }
 
 
 if __name__ == "__main__":
-    fileNames = ["cantilever_01", "portal_01", "grid_of_beams", "slab_01", "structure_01"]
+    fileNames = [
+        "cantilever_01",
+        "portal_01",
+        "grid_of_beams",
+        "slab_01",
+        "structure_01",
+        "building_02",
+    ]
     files = fileNames
 
     for fileName in files:

src/ifc2ca/scriptCodeAster.py

@@ -33,7 +33,9 @@ class COMMANDFILE:
             self.calculateRestraints(conn)
         for el in elements:
             for rel in el["connections"]:
-                conn = [c for c in connections if c["ifcName"] == rel["relatedConnection"]][0]
+                conn = [
+                    c for c in connections if c["ifcName"] == rel["relatedConnection"]
+                ][0]
                 rel["conn_string"] = None
                 if conn["geometryType"] == "point":
                     rel["conn_string"] = "_0DC_"
@@ -61,7 +63,10 @@ class COMMANDFILE:
                         + self.getGroupName(rel["relatingElement"])
                     )
                     rel["index"] = len(conn["relatedElements"]) + 1
-                    rel["unifiedGroupName"] = self.getGroupName(rel["relatedConnection"]) + "_0DC_%g" % rel["index"]
+                    rel["unifiedGroupName"] = (
+                        self.getGroupName(rel["relatedConnection"])
+                        + "_0DC_%g" % rel["index"]
+                    )
                 else:
                     rel["groupName2"] = self.getGroupName(rel["relatedConnection"])
                 self.calculateConstraints(rel)
@@ -71,25 +76,54 @@ class COMMANDFILE:
         materials = data["db"]["materials"]
         profiles = data["db"]["profiles"]
 
-        edgeGroupNames = tuple([self.getGroupName(el["ifcName"]) for el in elements if el["geometryType"] == "line"])
+        edgeGroupNames = tuple(
-        faceGroupNames = tuple([self.getGroupName(el["ifcName"]) for el in elements if el["geometryType"] == "surface"])
+            [
+                self.getGroupName(el["ifcName"])
+                for el in elements
+                if el["geometryType"] == "line"
+            ]
+        )
+        faceGroupNames = tuple(
+            [
+                self.getGroupName(el["ifcName"])
+                for el in elements
+                if el["geometryType"] == "surface"
+            ]
+        )
         point0DGroupNames = tuple(
-            [self.getGroupName(el["ifcName"]) + "_0D" for el in connections if el["geometryType"] == "point"]
+            [
+                self.getGroupName(el["ifcName"]) + "_0D"
+                for el in connections
+                if el["geometryType"] == "point"
+            ]
         )
         spring1DGroupNames = tuple(
             flatten(
-                [[rel["springGroupName"] for rel in el["connections"] if rel["springGroupName"]] for el in elements]
+                [
+                    [
+                        rel["springGroupName"]
+                        for rel in el["connections"]
+                        if rel["springGroupName"]
+                    ]
+                    for el in elements
+                ]
             )
         )
         point1DGroupNames = tuple(
-            [self.getGroupName(el["ifcName"]) + "_0D" for el in connections if el["geometryType"] == "line"]
+            [
+                self.getGroupName(el["ifcName"]) + "_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"]
+                rel["unifiedGroupName"]
+                for rel in conn["relatedElements"]
+                if rel["eccentricity"]
             ]
             # if not conn['appliedCondition'] and len(conn['unifiedGroupNames']) == 1:
             #     conn['appliedCondition'] = {
@@ -103,7 +137,9 @@ class COMMANDFILE:
                 unifiedConnection = True
             rigidLinkGroupNames.extend(
                 [
-                    self.getGroupName(rel["relatingElement"]) + "_1DR_" + self.getGroupName(conn["ifcName"])
+                    self.getGroupName(rel["relatingElement"])
+                    + "_1DR_"
+                    + self.getGroupName(conn["ifcName"])
                     for rel in conn["relatedElements"]
                     if rel["eccentricity"]
                 ]
@@ -182,7 +218,9 @@ model = AFFE_MODELE(
             MODELISATION = 'DIS_TR'
         ),"""
 
-            context = {"groupNames": tuple(flatten([point0DGroupNames, spring1DGroupNames]))}
+            context = {
+                "groupNames": tuple(flatten([point0DGroupNames, spring1DGroupNames]))
+            }
 
             f.write(template.format(**context))
 
@@ -234,7 +272,9 @@ model = AFFE_MODELE(
             else:
                 if "shearModulus" in material["mechProps"]:
                     poissonRatio = (
-                        material["mechProps"]["youngModulus"] / 2.0 / material["mechProps"]["shearModulus"]
+                        material["mechProps"]["youngModulus"]
+                        / 2.0
+                        / material["mechProps"]["shearModulus"]
                     ) - 1
                 else:
                     poissonRatio = 0.0
@@ -263,7 +303,9 @@ material = AFFE_MATERIAU(
         ),"""
 
             context = {
-                "groupNames": tuple([self.getGroupName(rel) for rel in material["relatedElements"]]),
+                "groupNames": tuple(
+                    [self.getGroupName(rel) for rel in material["relatedElements"]]
+                ),
                 "matNameID": "mat" + "_%s" % i,
             }
 
@@ -296,7 +338,10 @@ element = AFFE_CARA_ELEM(
         )
 
         for profile in profiles:
-            if profile["profileShape"] == "rectangular" and profile["profileType"] == "AREA":
+            if (
+                profile["profileShape"] == "rectangular"
+                and profile["profileType"] == "AREA"
+            ):
                 template = """
         _F(
             GROUP_MA = {groupNames},
@@ -306,13 +351,18 @@ element = AFFE_CARA_ELEM(
         ),"""
 
                 context = {
-                    "groupNames": tuple([self.getGroupName(rel) for rel in profile["relatedElements"]]),
+                    "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":
+            elif (
+                profile["profileShape"] == "iSymmetrical"
+                and profile["profileType"] == "AREA"
+            ):
                 template = """
         _F(
             GROUP_MA = {groupNames},
@@ -322,7 +372,9 @@ element = AFFE_CARA_ELEM(
         ),"""
 
                 context = {
-                    "groupNames": tuple([self.getGroupName(rel) for rel in profile["relatedElements"]]),
+                    "groupNames": tuple(
+                        [self.getGroupName(rel) for rel in profile["relatedElements"]]
+                    ),
                     "profileProperties": (
                         profile["mechProps"]["crossSectionArea"],
                         profile["mechProps"]["momentOfInertiaY"],
@@ -388,7 +440,10 @@ element = AFFE_CARA_ELEM(
             REPERE = 'LOCAL'
         ),"""
 
-            context = {"groupName": self.getGroupName(conn["ifcName"]) + "_0D", "stiffnesses": conn["stiffnesses"]}
+            context = {
+                "groupName": self.getGroupName(conn["ifcName"]) + "_0D",
+                "stiffnesses": conn["stiffnesses"],
+            }
 
             f.write(template.format(**context))
 
@@ -402,7 +457,10 @@ element = AFFE_CARA_ELEM(
             REPERE = 'LOCAL'
         ),"""
 
-                context = {"groupName": rel["springGroupName"], "stiffnesses": rel["stiffnesses"]}
+                context = {
+                    "groupName": rel["springGroupName"],
+                    "stiffnesses": rel["stiffnesses"],
+                }
 
                 f.write(template.format(**context))
 
@@ -416,7 +474,10 @@ element = AFFE_CARA_ELEM(
             REPERE = 'LOCAL'
         ),"""
 
-            context = {"groupName": self.getGroupName(conn["ifcName"]) + "_0D", "stiffnesses": conn["stiffnesses"]}
+            context = {
+                "groupName": self.getGroupName(conn["ifcName"]) + "_0D",
+                "stiffnesses": conn["stiffnesses"],
+            }
 
             f.write(template.format(**context))
 
@@ -439,7 +500,10 @@ element = AFFE_CARA_ELEM(
             VALE = {localAxisY}
         ),"""
 
-            context = {"groupName": self.getGroupName(el["ifcName"]), "localAxisY": tuple(el["orientation"][1])}
+            context = {
+                "groupName": self.getGroupName(el["ifcName"]),
+                "localAxisY": tuple(el["orientation"][1]),
+            }
 
             f.write(template.format(**context))
 
@@ -614,7 +678,9 @@ liaisons = AFFE_CHAR_MECA(
     LIAISON_UNIF = ("""
             )
 
-            for conn in [conn for conn in connections if len(conn["unifiedGroupNames"]) > 1]:
+            for conn in [
+                conn for conn in connections if len(conn["unifiedGroupNames"]) > 1
+            ]:
                 template = """
         _F(
             GROUP_NO = {groupNames},
@@ -798,44 +864,103 @@ FIN()
         gr1 = rel["groupName1"]
         gr2 = rel["groupName2"]
         o = np.array(rel["orientation"]).transpose().tolist()
-        liaisons = {"groupNames": (gr1, gr1, gr1, gr2, gr2, gr2), "coeffs": [], "dofs": []}
+        liaisons = {
+            "groupNames": (gr1, gr1, gr1, gr2, gr2, gr2),
+            "coeffs": [],
+            "dofs": [],
+        }
         stiffnesses = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
         if not rel["appliedCondition"]:
-            rel["appliedCondition"] = {"dx": True, "dy": True, "dz": True, "drx": True, "dry": True, "drz": True}
+            rel["appliedCondition"] = {
-        if isinstance(rel["appliedCondition"]["dx"], bool) and rel["appliedCondition"]["dx"]:
+                "dx": True,
-            liaisons["coeffs"].append((o[0][0], o[1][0], o[2][0], -o[0][0], -o[1][0], -o[2][0]))
+                "dy": True,
+                "dz": True,
+                "drx": True,
+                "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])
+            )
             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"]:
+        if (
-            liaisons["coeffs"].append((o[0][1], o[1][1], o[2][1], -o[0][1], -o[1][1], -o[2][1]))
+            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"]:
+        if (
-            liaisons["coeffs"].append((o[0][2], o[1][2], o[2][2], -o[0][2], -o[1][2], -o[2][2]))
+            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"]:
+        if (
-            liaisons["coeffs"].append((o[0][0], o[1][0], o[2][0], -o[0][0], -o[1][0], -o[2][0]))
+            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"]:
+        if (
-            liaisons["coeffs"].append((o[0][1], o[1][1], o[2][1], -o[0][1], -o[1][1], -o[2][1]))
+            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"]:
+        if (
-            liaisons["coeffs"].append((o[0][2], o[1][2], o[2][2], -o[0][2], -o[1][2], -o[2][2]))
+            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
@@ -852,40 +977,76 @@ 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
@@ -893,7 +1054,13 @@ FIN()
 
 
 if __name__ == "__main__":
-    fileNames = ["cantilever_01", "portal_01", "grid_of_beams", "slab_01", "structure_01"]
+    fileNames = [
+        "cantilever_01",
+        "portal_01",
+        "grid_of_beams",
+        "slab_01",
+        "structure_01",
+    ]
     files = fileNames
 
     for fileName in files:

src/ifc2ca/scriptCodeAsterBonded.py

@@ -0,0 +1,555 @@
+import json
+import numpy as np
+import itertools
+
+flatten = itertools.chain.from_iterable
+
+ScaleFactor = 1.0
+
+AccelOfGravity = 9.806 * 1000
+
+class COMMANDFILE:
+    def __init__(self, dataFilename, asterFilename):
+        self.dataFilename = dataFilename
+        self.asterFilename = asterFilename
+        self.create()
+
+    def getGroupName(self, name):
+        info = name.split("|")
+        sortName = "".join(c for c in info[0] if c.isupper())
+        return str(sortName + "_" + info[1])
+
+    def create(self):
+        # 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"] = []
+        for el in elements:
+            for rel in el["connections"]:
+                conn = [
+                    c for c in connections if c["ifcName"] == rel["relatedConnection"]
+                ][0]
+                conn["relatedElements"].append(rel)
+        # End <--
+
+        materials = data["db"]["materials"]
+        profiles = data["db"]["profiles"]
+
+        edgeGroupNames = tuple(
+            [
+                self.getGroupName(el["ifcName"])
+                for el in elements
+                if el["geometryType"] == "line"
+            ]
+        )
+        faceGroupNames = tuple(
+            [
+                self.getGroupName(el["ifcName"])
+                for el in elements
+                if el["geometryType"] == "surface"
+            ]
+        )
+
+        rigidLinkGroupNames = []
+        for conn in connections:
+            rigidLinkGroupNames.extend(
+                [
+                    self.getGroupName(rel["relatingElement"])
+                    + "_1DR_"
+                    + self.getGroupName(conn["ifcName"])
+                    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 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"])
+                * ScaleFactor ** 2,
+                "poissonRatio": float(poissonRatio),
+                "massDensity": float(material["commonProps"]["massDensity"])
+                * ScaleFactor ** 3,
+            }
+
+            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"] / ScaleFactor,
+                        profile["yDim"] / ScaleFactor,
+                    ),
+                }
+
+                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"] / ScaleFactor ** 2,
+                        profile["mechProps"]["momentOfInertiaY"] / ScaleFactor ** 4,
+                        profile["mechProps"]["momentOfInertiaZ"] / ScaleFactor ** 4,
+                        profile["mechProps"]["torsionalConstantX"] / ScaleFactor ** 4,
+                    ),
+                }
+
+                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["ifcName"]),
+                "thickness": el["thickness"] / ScaleFactor,
+                "localAxisX": tuple(el["orientation"][0]),
+            }
+
+            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["ifcName"]),
+                "localAxisY": tuple(el["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,
+    DDL_IMPO = (
+        _F(
+            GROUP_NO = 'grdSupps',
+            DX = 0.0,
+            DY = 0.0,
+            DZ = 0.0,
+            DRX = 0.0,
+            DRY = 0.0,
+            DRZ = 0.0
+        )
+    ),"""
+        )
+
+        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()
+
+
+if __name__ == "__main__":
+    fileNames = ["building_02"]
+    files = fileNames
+
+    for fileName in files:
+        BASE_PATH = "/home/jesusbill/Dev-Projects/github.com/IfcOpenShell/analysis-models/models/"
+        DATAFILENAME = BASE_PATH + fileName + "/" + fileName + ".json"
+        ASTERFILENAME = BASE_PATH + fileName + "/" + fileName + ".comm"
+        COMMANDFILE(DATAFILENAME, ASTERFILENAME)

src/ifc2ca/scriptSalome.py

@@ -75,7 +75,9 @@ class MODEL:
             shapeType = "EDGE"
         if geometryType == "surface":
             shapeType = "FACE"
-        return self.geompy.MakePartition(objects, [], [], [], self.geompy.ShapeType[shapeType], 0, [], 1)
+        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"])
@@ -157,17 +159,25 @@ class MODEL:
             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["ifcName"] == rel["relatedConnection"]][0]
+                conn = [
+                    c for c in connections if c["ifcName"] == rel["relatedConnection"]
+                ][0]
                 if rel["eccentricity"]:
                     rel["index"] = len(conn["relatedElements"]) + 1
                 conn["relatedElements"].append(rel)
 
                 if not rel["eccentricity"]:
-                    el["connObjs"][j] = self.makeObject(conn["geometry"], conn["geometryType"])
+                    el["connObjs"][j] = self.makeObject(
+                        conn["geometry"], conn["geometryType"]
+                    )
                 else:
                     if conn["geometryType"] == "point":
-                        geometry = self.getLinkGeometry(rel["eccentricity"], el["orientation"], conn["geometry"])
+                        geometry = self.getLinkGeometry(
-                        el["connObjs"][j] = self.makeObject(geometry[0], conn["geometryType"])
+                            rel["eccentricity"], el["orientation"], conn["geometry"]
+                        )
+                        el["connObjs"][j] = self.makeObject(
+                            geometry[0], conn["geometryType"]
+                        )
 
                         el["linkPointObjs"][j][0] = self.geompy.MakeVertex(
                             geometry[0][0], geometry[0][1], geometry[0][2]
@@ -179,9 +189,14 @@ class MODEL:
                             el["linkPointObjs"][j][0], el["linkPointObjs"][j][1]
                         )
                     else:
-                        print("Eccentricity defined for a %s geometryType" % conn["geometryType"])
+                        print(
+                            "Eccentricity defined for a %s geometryType"
+                            % conn["geometryType"]
+                        )
 
-            el["partObj"] = self.makePartition([el["elemObj"]] + el["connObjs"], el["geometryType"])
+            el["partObj"] = self.makePartition(
+                [el["elemObj"]] + el["connObjs"], el["geometryType"]
+            )
 
             el["elemObj"] = geompy.GetInPlace(el["partObj"], el["elemObj"])
             for j, rel in enumerate(el["connections"]):
@@ -193,7 +208,9 @@ class MODEL:
         # 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)
@@ -203,9 +220,13 @@ class MODEL:
         # Loop 2
         for el in elements:
             # geompy.addToStudy(el['partObj'], self.getGroupName(el['ifcName']))
-            geompy.addToStudyInFather(el["partObj"], el["elemObj"], self.getGroupName(el["ifcName"]))
+            geompy.addToStudyInFather(
+                el["partObj"], el["elemObj"], self.getGroupName(el["ifcName"])
+            )
             for j, rel in enumerate(el["connections"]):
-                conn = [c for c in connections if c["ifcName"] == rel["relatedConnection"]][0]
+                conn = [
+                    c for c in connections if c["ifcName"] == rel["relatedConnection"]
+                ][0]
                 rel["conn_string"] = None
                 if conn["geometryType"] == "point":
                     rel["conn_string"] = "_0DC_"
@@ -216,7 +237,9 @@ class MODEL:
                 geompy.addToStudyInFather(
                     el["partObj"],
                     el["connObjs"][j],
-                    self.getGroupName(el["ifcName"]) + rel["conn_string"] + self.getGroupName(rel["relatedConnection"]),
+                    self.getGroupName(el["ifcName"])
+                    + rel["conn_string"]
+                    + self.getGroupName(rel["relatedConnection"]),
                 )
                 if rel["eccentricity"]:
                     pass
@@ -226,7 +249,9 @@ class MODEL:
 
         for conn in connections:
             # geompy.addToStudy(conn['connObj'], self.getGroupName(conn['ifcName']))
-            geompy.addToStudyInFather(conn["connObj"], conn["connObj"], self.getGroupName(conn["ifcName"]))
+            geompy.addToStudyInFather(
+                conn["connObj"], conn["connObj"], self.getGroupName(conn["ifcName"])
+            )
 
         elapsed_time = time.time() - init_time
         init_time += elapsed_time
@@ -240,14 +265,18 @@ class MODEL:
         # Define and add groups for all curve and surface members
         if len([e for e in elements if e["geometryType"] == "line"]) > 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["geometryType"] == "line"]
+            )
             # Define group object and add to study
             curveCompound = geompy.GetInPlace(bldComp, compoundTemp)
             geompy.addToStudyInFather(bldComp, curveCompound, "CurveMembers")
 
         if len([e for e in elements if e["geometryType"] == "surface"]) > 0:
             # Make compound of requested group
-            compoundTemp = geompy.MakeCompound([e["elemObj"] for e in elements if e["geometryType"] == "surface"])
+            compoundTemp = geompy.MakeCompound(
+                [e["elemObj"] for e in elements if e["geometryType"] == "surface"]
+            )
             # Define group object and add to study
             surfaceCompound = geompy.GetInPlace(bldComp, compoundTemp)
             geompy.addToStudyInFather(bldComp, surfaceCompound, "SurfaceMembers")
@@ -255,34 +284,45 @@ 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["ifcName"]))
+            geompy.addToStudyInFather(
+                bldComp, el["elemObj"], self.getGroupName(el["ifcName"])
+            )
 
             for j, rel in enumerate(el["connections"]):
                 geompy.addToStudyInFather(
                     bldComp,
                     el["connObjs"][j],
-                    self.getGroupName(el["ifcName"]) + rel["conn_string"] + self.getGroupName(rel["relatedConnection"]),
+                    self.getGroupName(el["ifcName"])
+                    + rel["conn_string"]
+                    + self.getGroupName(rel["relatedConnection"]),
                 )
                 if rel["eccentricity"]:  # point geometry
                     geompy.addToStudyInFather(
                         bldComp,
                         el["linkObjs"][j],
-                        self.getGroupName(el["ifcName"]) + "_1DR_" + self.getGroupName(rel["relatedConnection"]),
+                        self.getGroupName(el["ifcName"])
+                        + "_1DR_"
+                        + self.getGroupName(rel["relatedConnection"]),
                     )
                     geompy.addToStudyInFather(
                         bldComp,
                         el["linkPointObjs"][j][0],
-                        self.getGroupName(rel["relatedConnection"]) + "_0DC_" + self.getGroupName(el["ifcName"]),
+                        self.getGroupName(rel["relatedConnection"])
+                        + "_0DC_"
+                        + self.getGroupName(el["ifcName"]),
                     )
                     geompy.addToStudyInFather(
                         bldComp,
                         el["linkPointObjs"][j][1],
-                        self.getGroupName(rel["relatedConnection"]) + "_0DC_%g" % rel["index"],
+                        self.getGroupName(rel["relatedConnection"])
+                        + "_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["ifcName"]))
+            geompy.addToStudyInFather(
+                bldComp, conn["connObj"], self.getGroupName(conn["ifcName"])
+            )
 
         elapsed_time = time.time() - init_time
         init_time += elapsed_time
@@ -339,7 +379,9 @@ class MODEL:
             smesh.SetName(tempgroup, "CurveMembers")
 
         if len([e for e in elements if e["geometryType"] == "surface"]) > 0:
-            tempgroup = bldMesh.GroupOnGeom(surfaceCompound, "SurfaceMembers", SMESH.FACE)
+            tempgroup = bldMesh.GroupOnGeom(
+                surfaceCompound, "SurfaceMembers", SMESH.FACE
+            )
             smesh.SetName(tempgroup, "SurfaceMembers")
 
         # Define groups in Mesh
@@ -348,41 +390,59 @@ class MODEL:
                 shapeType = SMESH.EDGE
             if el["geometryType"] == "surface":
                 shapeType = SMESH.FACE
-            tempgroup = bldMesh.GroupOnGeom(el["elemObj"], self.getGroupName(el["ifcName"]), shapeType)
+            tempgroup = bldMesh.GroupOnGeom(
+                el["elemObj"], self.getGroupName(el["ifcName"]), shapeType
+            )
             smesh.SetName(tempgroup, self.getGroupName(el["ifcName"]))
 
             for j, rel in enumerate(el["connections"]):
                 tempgroup = bldMesh.GroupOnGeom(
                     el["connObjs"][j],
-                    self.getGroupName(el["ifcName"]) + rel["conn_string"] + self.getGroupName(rel["relatedConnection"]),
+                    self.getGroupName(el["ifcName"])
+                    + rel["conn_string"]
+                    + self.getGroupName(rel["relatedConnection"]),
                     SMESH.NODE,
                 )
                 smesh.SetName(
                     tempgroup,
-                    self.getGroupName(el["ifcName"]) + rel["conn_string"] + self.getGroupName(rel["relatedConnection"]),
+                    self.getGroupName(el["ifcName"])
+                    + rel["conn_string"]
+                    + self.getGroupName(rel["relatedConnection"]),
                 )
-                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["ifcName"]) + "_1DR_" + self.getGroupName(rel["relatedConnection"]),
+                        self.getGroupName(el["ifcName"])
+                        + "_1DR_"
+                        + self.getGroupName(rel["relatedConnection"]),
                         SMESH.EDGE,
                     )
                     smesh.SetName(
                         tempgroup,
-                        self.getGroupName(el["ifcName"]) + "_1DR_" + self.getGroupName(rel["relatedConnection"]),
+                        self.getGroupName(el["ifcName"])
+                        + "_1DR_"
+                        + self.getGroupName(rel["relatedConnection"]),
                     )
 
                     tempgroup = bldMesh.GroupOnGeom(
                         el["linkPointObjs"][j][0],
-                        self.getGroupName(rel["relatedConnection"]) + "_0DC_" + self.getGroupName(el["ifcName"]),
+                        self.getGroupName(rel["relatedConnection"])
+                        + "_0DC_"
+                        + self.getGroupName(el["ifcName"]),
                         SMESH.NODE,
                     )
                     smesh.SetName(
                         tempgroup,
-                        self.getGroupName(rel["relatedConnection"]) + "_0DC_" + self.getGroupName(el["ifcName"]),
+                        self.getGroupName(rel["relatedConnection"])
+                        + "_0DC_"
+                        + self.getGroupName(el["ifcName"]),
                     )
-                    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],
@@ -391,38 +451,60 @@ class MODEL:
                         + self.getGroupName(rel["relatedConnection"]),
                         SMESH.NODE,
                     )
-                    smesh.SetName(tempgroup, self.getGroupName(rel["relatedConnection"]) + "_0DC_%g" % rel["index"])
+                    smesh.SetName(
+                        tempgroup,
+                        self.getGroupName(rel["relatedConnection"])
+                        + "_0DC_%g" % rel["index"],
+                    )
 
         for conn in connections:
-            tempgroup = bldMesh.GroupOnGeom(conn["connObj"], self.getGroupName(conn["ifcName"]), SMESH.NODE)
+            tempgroup = bldMesh.GroupOnGeom(
+                conn["connObj"], self.getGroupName(conn["ifcName"]), SMESH.NODE
+            )
             smesh.SetName(tempgroup, self.getGroupName(conn["ifcName"]))
             nodesId = bldMesh.GetIDSource(tempgroup.GetNodeIDs(), SMESH.NODE)
-            tempgroup = bldMesh.Add0DElementsToAllNodes(nodesId, self.getGroupName(conn["ifcName"]))
+            tempgroup = bldMesh.Add0DElementsToAllNodes(
+                nodesId, self.getGroupName(conn["ifcName"])
+            )
             smesh.SetName(tempgroup, self.getGroupName(conn["ifcName"] + "_0D"))
             if conn["geometryType"] == "point":
                 conn["node"] = nodesId.GetIDs()[0]
             if conn["geometryType"] == "line":
-                tempgroup = bldMesh.GroupOnGeom(conn["connObj"], self.getGroupName(conn["ifcName"]), SMESH.EDGE)
+                tempgroup = bldMesh.GroupOnGeom(
+                    conn["connObj"], self.getGroupName(conn["ifcName"]), SMESH.EDGE
+                )
                 smesh.SetName(tempgroup, self.getGroupName(conn["ifcName"]))
             if conn["geometryType"] == "surface":
-                tempgroup = bldMesh.GroupOnGeom(conn["connObj"], self.getGroupName(conn["ifcName"]), SMESH.FACE)
+                tempgroup = bldMesh.GroupOnGeom(
+                    conn["connObj"], self.getGroupName(conn["ifcName"]), SMESH.FACE
+                )
                 smesh.SetName(tempgroup, self.getGroupName(conn["ifcName"]))
 
         # create 1D SEG2 spring elements
         for el in elements:
             for j, rel in enumerate(el["connections"]):
-                conn = [c for c in connections if c["ifcName"] == rel["relatedConnection"]][0]
+                conn = [
+                    c for c in connections if c["ifcName"] == rel["relatedConnection"]
+                ][0]
                 if conn["geometryType"] == "point":
                     grpName = bldMesh.CreateEmptyGroup(
                         SMESH.EDGE,
-                        self.getGroupName(el["ifcName"]) + "_1DS_" + self.getGroupName(rel["relatedConnection"]),
+                        self.getGroupName(el["ifcName"])
+                        + "_1DS_"
+                        + self.getGroupName(rel["relatedConnection"]),
                     )
                     smesh.SetName(
                         grpName,
-                        self.getGroupName(el["ifcName"]) + "_1DS_" + self.getGroupName(rel["relatedConnection"]),
+                        self.getGroupName(el["ifcName"])
+                        + "_1DS_"
+                        + self.getGroupName(rel["relatedConnection"]),
                     )
                     if not rel["eccentricity"]:
-                        conn = [conn for conn in connections if conn["ifcName"] == rel["relatedConnection"]][0]
+                        conn = [
+                            conn
+                            for conn in connections
+                            if conn["ifcName"] == rel["relatedConnection"]
+                        ][0]
                         grpName.Add([bldMesh.AddEdge([conn["node"], rel["node"]])])
                     else:
                         grpName.Add([bldMesh.AddEdge([rel["eccNode"], rel["node"]])])
@@ -437,7 +519,12 @@ class MODEL:
         try:
             if NEW_SALOME:
                 bldMesh.ExportMED(
-                    self.medFilename, auto_groups=0, minor=40, overwrite=1, meshPart=None, autoDimension=0
+                    self.medFilename,
+                    auto_groups=0,
+                    minor=40,
+                    overwrite=1,
+                    meshPart=None,
+                    autoDimension=0,
                 )
             else:
                 bldMesh.ExportMED(self.medFilename, 0, SMESH.MED_V2_2, 1, None, 0)

src/ifc2ca/scriptSalomeBonded.py

@@ -0,0 +1,415 @@
+from __future__ import division
+from __future__ import print_function
+import os
+import time
+import json
+import salome
+import salome_notebook
+import salome_version
+import numpy as np
+import itertools
+
+flatten = itertools.chain.from_iterable
+
+ScaleFactor = 1.0
+decimals = 2
+
+
+class MODEL:
+    def __init__(self, dataFilename, medFilename, meshSize, zGround):
+        self.dataFilename = dataFilename
+        self.medFilename = medFilename
+        self.meshSize = meshSize
+        self.zGround = zGround
+        self.tolLoc = 0
+        self.mesh = None
+        self.meshNodes = None
+        self.create()
+
+    def getGroupName(self, name):
+        info = name.split("|")
+        sortName = "".join(c for c in info[0] if c.isupper())
+        return str(sortName + "_" + info[1])
+
+    def makePoint(self, pl):
+        """Function to define a Point from
+        a polyline (list of 1 point)"""
+
+        (x, y, z) = pl
+        return self.geompy.MakeVertex(x, y, z)
+
+    def makeLine(self, pl):
+        """Function to define a Line from
+        a polyline (list of 2 points)"""
+
+        (x, y, z) = pl[0]
+        P1 = self.geompy.MakeVertex(x, y, z)
+        (x, y, z) = pl[1]
+        P2 = self.geompy.MakeVertex(x, y, z)
+
+        return self.geompy.MakeLineTwoPnt(P1, P2)
+
+    def makeFace(self, pl):
+        """Function to define a Face from
+        a polyline (list of points)"""
+
+        pointList = [None for _ in range(len(pl))]
+        for ip, (x, y, z) in enumerate(pl):
+            pointList[ip] = self.geompy.MakeVertex(x, y, z)
+
+        LineList = [None for _ in range(len(pl))]
+        for ip, P2 in enumerate(pointList):
+            P1 = pointList[ip - 1]
+            LineList[ip] = self.geompy.MakeLineTwoPnt(P1, P2)
+
+        return self.geompy.MakeFaceWires(LineList, 1)
+
+    def makeObject(self, geometry, geometryType):
+        geometry = np.round(np.array(geometry), decimals=decimals) / ScaleFactor
+        geometry = geometry.tolist()
+        if geometryType == "point":
+            return self.makePoint(geometry)
+        if geometryType == "line":
+            return self.makeLine(geometry)
+        if geometryType == "surface":
+            return self.makeFace(geometry)
+
+    def makePartition(self, objects, geometryType):
+        if geometryType == "point":
+            shapeType = "VERTEX"
+        if geometryType == "line":
+            shapeType = "EDGE"
+        if geometryType == "surface":
+            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]
+
+    def length(self, geometry):
+        return (
+            (geometry[1][0] - geometry[0][0]) ** 2
+            + (geometry[1][1] - geometry[0][1]) ** 2
+            + (geometry[1][2] - geometry[0][2]) ** 2
+        ) ** 0.5
+
+    def select(self, elements):
+        zmin = -100
+        zmax = 126300
+        for el in elements:
+            include = True
+            for p in el["geometry"]:
+                if p[2] < zmin or p[2] > zmax:
+                    include = False
+                    break
+            el["include"] = include
+        return [el for el in elements if el["include"]]
+
+    def create(self):
+        # Read data from input file
+        with open(self.dataFilename) as dataFile:
+            data = json.load(dataFile)
+
+        # print(len(data['elements']))
+        # elements = self.select(data['elements'])
+        # print(len(elements))
+        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"] = []
+        # End <--
+
+        meshSize = self.meshSize / ScaleFactor
+        zGround = self.zGround / ScaleFactor
+
+        dec = 5  # 4 decimals for length in mm
+        tol = 10 ** (-dec - 3 + 1)
+
+        self.tolLoc = tol * 10 * 2
+        tolLoc = self.tolLoc
+
+        NEW_SALOME = int(salome_version.getVersion()[0]) >= 9
+        salome.salome_init()
+        theStudy = salome.myStudy
+        notebook = salome_notebook.NoteBook(theStudy)
+
+        ###
+        ### GEOM component
+        ###
+        import GEOM
+        from salome.geom import geomBuilder
+        import math
+        import SALOMEDS
+
+        gg = salome.ImportComponentGUI("GEOM")
+        if NEW_SALOME:
+            geompy = geomBuilder.New()
+        else:
+            geompy = geomBuilder.New(theStudy)
+        self.geompy = geompy
+
+        O = geompy.MakeVertex(0, 0, 0)
+        OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
+        OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
+        OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)
+        geompy.addToStudy(O, "O")
+        geompy.addToStudy(OX, "OX")
+        geompy.addToStudy(OY, "OY")
+        geompy.addToStudy(OZ, "OZ")
+
+        if len([e for e in elements if e["geometryType"] == "line"]) > 0:
+            buildingShapeType = "EDGE"
+        if len([e for e in elements if e["geometryType"] == "surface"]) > 0:
+            buildingShapeType = "FACE"
+
+        ### Define entities ###
+        start_time = time.time()
+        print("Defining Object Geometry")
+        init_time = start_time
+
+        # Loop 1
+        for el in elements:
+            el["elemObj"] = self.makeObject(el["geometry"], el["geometryType"])
+
+            el["linkObjs"] = [None for _ in el["connections"]]
+            for j, rel in enumerate(el["connections"]):
+                conn = [
+                    c for c in connections if c["ifcName"] == rel["relatedConnection"]
+                ][0]
+                if rel["eccentricity"]:
+                    rel["index"] = len(conn["relatedElements"]) + 1
+
+                    geometry = self.getLinkGeometry(
+                        rel["eccentricity"], el["orientation"], conn["geometry"]
+                    )
+                    el["linkObjs"][j] = self.makeObject(geometry, "line")
+                conn["relatedElements"].append(rel)
+
+        # Make assemble of Building Object
+        bldObjs = []
+        bldObjs.extend([el["elemObj"] for el in elements])
+        bldObjs.extend(
+            flatten([[link for link in el["linkObjs"] if link] for el in elements])
+        )
+
+        # bldComp = geompy.MakeCompound(bldObjs)
+        bldComp = geompy.MakePartition(
+            bldObjs, [], [], [], self.geompy.ShapeType[buildingShapeType], 0, [], 1
+        )
+        geompy.addToStudy(bldComp, "bldComp")
+
+        elapsed_time = time.time() - init_time
+        init_time += elapsed_time
+        print("Building Geometry Defined in %g sec" % (elapsed_time))
+
+        if len([e for e in elements if e["geometryType"] == "line"]) > 0:
+            buildingShapeType = "EDGE"
+        if len([e for e in elements if e["geometryType"] == "surface"]) > 0:
+            buildingShapeType = "FACE"
+
+        # Define and add groups for all curve, surface and rigid members
+        if len([e for e in elements if e["geometryType"] == "line"]) > 0:
+            # Make compound of requested group
+            compoundTemp = geompy.MakeCompound(
+                [e["elemObj"] for e in elements if e["geometryType"] == "line"]
+            )
+            # Define group object and add to study
+            curveCompound = geompy.GetInPlace(bldComp, compoundTemp)
+            geompy.addToStudyInFather(bldComp, curveCompound, "CurveMembers")
+
+        if len([e for e in elements if e["geometryType"] == "surface"]) > 0:
+            # Make compound of requested group
+            compoundTemp = geompy.MakeCompound(
+                [e["elemObj"] for e in elements if e["geometryType"] == "surface"]
+            )
+            # Define group object and add to study
+            surfaceCompound = geompy.GetInPlace(bldComp, compoundTemp)
+            geompy.addToStudyInFather(bldComp, surfaceCompound, "SurfaceMembers")
+
+        linkObjs = list(
+            flatten([[obj for obj in el["linkObjs"] if obj] for el in elements])
+        )
+        if len(linkObjs) > 0:
+            # Make compound of requested group
+            compoundTemp = geompy.MakeCompound(linkObjs)
+            # Define group object and add to study
+            rigidCompound = geompy.GetInPlace(bldComp, compoundTemp)
+            geompy.addToStudyInFather(bldComp, rigidCompound, "RigidMembers")
+
+        for el in elements:
+            # el['partObj'] = geompy.RestoreGivenSubShapes(bldComp, [el['partObj']], GEOM.FSM_GetInPlace, False, False)[0]
+            el["elemObj"] = geompy.GetInPlace(bldComp, el["elemObj"])
+            geompy.addToStudyInFather(
+                bldComp, el["elemObj"], self.getGroupName(el["ifcName"])
+            )
+
+            for j, rel in enumerate(el["connections"]):
+                if rel["eccentricity"]:  # point geometry
+                    el["linkObjs"][j] = geompy.GetInPlace(bldComp, el["linkObjs"][j])
+                    geompy.addToStudyInFather(
+                        bldComp,
+                        el["linkObjs"][j],
+                        self.getGroupName(el["ifcName"])
+                        + "_1DR_"
+                        + self.getGroupName(rel["relatedConnection"]),
+                    )
+
+        elapsed_time = time.time() - init_time
+        init_time += elapsed_time
+        print("Building Geometry Groups Defined in %g sec" % (elapsed_time))
+
+        ###
+        ### SMESH component
+        ###
+
+        import SMESH
+        from salome.smesh import smeshBuilder
+
+        print("Defining Mesh Components")
+
+        if NEW_SALOME:
+            smesh = smeshBuilder.New()
+        else:
+            smesh = smeshBuilder.New(theStudy)
+        bldMesh = smesh.Mesh(bldComp)
+        Regular_1D = bldMesh.Segment()
+        Local_Length_1 = Regular_1D.LocalLength(meshSize, None, tolLoc)
+
+        if buildingShapeType == "FACE":
+            NETGEN2D_ONLY = bldMesh.Triangle(algo=smeshBuilder.NETGEN_2D)
+            NETGEN2D_Pars = NETGEN2D_ONLY.Parameters()
+            NETGEN2D_Pars.SetMaxSize(meshSize)
+            NETGEN2D_Pars.SetOptimize(1)
+            NETGEN2D_Pars.SetFineness(2)
+            NETGEN2D_Pars.SetMinSize(meshSize / 5.0)
+            NETGEN2D_Pars.SetUseSurfaceCurvature(1)
+            NETGEN2D_Pars.SetQuadAllowed(1)
+            NETGEN2D_Pars.SetSecondOrder(0)
+            NETGEN2D_Pars.SetFuseEdges(254)
+
+        isDone = bldMesh.Compute()
+        coincident_nodes_on_part = bldMesh.FindCoincidentNodesOnPart(
+            [bldMesh], tolLoc, [], 0
+        )
+        if coincident_nodes_on_part:
+            # bldMesh.MergeNodes(coincident_nodes_on_part, [], 0)
+            # print(f'{len(coincident_nodes_on_part)} Sets of Coincident Nodes Found and Merged')
+            print(f"{len(coincident_nodes_on_part)} Sets of Coincident Nodes Found")
+            print(f"{coincident_nodes_on_part}")
+
+        ## Set names of Mesh objects
+        smesh.SetName(Regular_1D.GetAlgorithm(), "Regular_1D")
+        smesh.SetName(Local_Length_1, "Local_Length_1")
+
+        if buildingShapeType == "FACE":
+            smesh.SetName(NETGEN2D_ONLY.GetAlgorithm(), "NETGEN2D_ONLY")
+            smesh.SetName(NETGEN2D_Pars, "NETGEN2D_Pars")
+
+        smesh.SetName(bldMesh.GetMesh(), "bldMesh")
+
+        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, surface and rigid members
+        if len([e for e in elements if e["geometryType"] == "line"]) > 0:
+            tempgroup = bldMesh.GroupOnGeom(curveCompound, "CurveMembers", SMESH.EDGE)
+            smesh.SetName(tempgroup, "CurveMembers")
+
+        if len([e for e in elements if e["geometryType"] == "surface"]) > 0:
+            tempgroup = bldMesh.GroupOnGeom(
+                surfaceCompound, "SurfaceMembers", SMESH.FACE
+            )
+            smesh.SetName(tempgroup, "SurfaceMembers")
+
+        if len(linkObjs) > 0:
+            tempgroup = bldMesh.GroupOnGeom(rigidCompound, "RigidMembers", SMESH.EDGE)
+            smesh.SetName(tempgroup, "RigidMembers")
+
+        # Define groups in Mesh
+        for el in elements:
+            if el["geometryType"] == "line":
+                shapeType = SMESH.EDGE
+            if el["geometryType"] == "surface":
+                shapeType = SMESH.FACE
+            tempgroup = bldMesh.GroupOnGeom(
+                el["elemObj"], self.getGroupName(el["ifcName"]), shapeType
+            )
+            smesh.SetName(tempgroup, self.getGroupName(el["ifcName"]))
+
+            for j, rel in enumerate(el["connections"]):
+                if rel["eccentricity"]:
+                    tempgroup = bldMesh.GroupOnGeom(
+                        el["linkObjs"][j],
+                        self.getGroupName(el["ifcName"])
+                        + "_1DR_"
+                        + self.getGroupName(rel["relatedConnection"]),
+                        SMESH.EDGE,
+                    )
+                    smesh.SetName(
+                        tempgroup,
+                        self.getGroupName(el["ifcName"])
+                        + "_1DR_"
+                        + self.getGroupName(rel["relatedConnection"]),
+                    )
+
+        self.mesh = bldMesh
+        self.meshNodes = bldMesh.GetNodesId()
+
+        # Find ground supports and extract node coordinates
+        grdSupps = bldMesh.CreateEmptyGroup(SMESH.NODE, "grdSupps")
+
+        for node in self.meshNodes:
+            coords = bldMesh.GetNodeXYZ(node)
+            if abs(coords[2] - self.zGround) < tolLoc:
+                grdSupps.Add([node])
+
+        smesh.SetName(grdSupps, "grdSupps")
+
+        elapsed_time = time.time() - init_time
+        init_time += elapsed_time
+        print("Mesh Groups Defined in %g sec" % (elapsed_time))
+
+        try:
+            if NEW_SALOME:
+                bldMesh.ExportMED(
+                    self.medFilename,
+                    auto_groups=0,
+                    minor=40,
+                    overwrite=1,
+                    meshPart=None,
+                    autoDimension=0,
+                )
+            else:
+                bldMesh.ExportMED(self.medFilename, 0, SMESH.MED_V2_2, 1, None, 0)
+        except:
+            print("ExportMED() failed. Invalid file name?")
+
+        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 = ["building_02"]
+    files = fileNames
+
+    meshSize = 500
+    zGround = 0
+
+    for fileName in files:
+        BASE_PATH = "/home/jesusbill/Dev-Projects/github.com/IfcOpenShell/analysis-models/models/"
+        DATAFILENAME = BASE_PATH + fileName + "/" + fileName + ".json"
+        MEDFILENAME = BASE_PATH + fileName + "/" + fileName + ".med"
+        model = MODEL(DATAFILENAME, MEDFILENAME, meshSize, zGround)