Force creation of point elements

kratos.gid/kratos.tcl

@@ -9,7 +9,6 @@ namespace eval ::Kratos {
     variable must_write_calc_data
     variable must_exist_calc_data
 
-    variable tmp_init_mesh_time
     variable namespaces
 
     variable mesh_criteria_forced
@@ -200,7 +199,7 @@ proc Kratos::LoadCommonScripts { } {
         uplevel #0 [list source [file join $kratos_private(Path) scripts Writing $filename.tcl]]
     }
     # Common scripts
-    foreach filename {Utils Launch Applications spdAuxiliar Menus Deprecated Logs} {
+    foreach filename {Utils Launch Applications spdAuxiliar Mesh Menus Deprecated Logs} {
         uplevel #0 [list source [file join $kratos_private(Path) scripts $filename.tcl]]
     }
     # Common controllers
@@ -427,19 +426,10 @@ proc Kratos::Event_BeforeMeshGeneration {elementsize} {
     set tmp_init_mesh_time $inittime
     Kratos::Log "Mesh BeforeMeshGeneration start"
 
-    GiD_MeshData mesh_criteria to_be_meshed 1 lines [GiD_Geometry list line]
-    GiD_MeshData mesh_criteria to_be_meshed 1 surfaces [GiD_Geometry list surface]
-    GiD_MeshData mesh_criteria to_be_meshed 1 volumes  [GiD_Geometry list volume ]
-
-    # We need to mesh every line and surface assigned to a group that appears in the tree
-    foreach group [spdAux::GetAppliedGroups] {
-        GiD_MeshData mesh_criteria to_be_meshed 2 lines [GiD_EntitiesGroups get $group lines]
-        GiD_MeshData mesh_criteria to_be_meshed 2 surfaces [GiD_EntitiesGroups get $group surfaces]
-        GiD_MeshData mesh_criteria to_be_meshed 2 volumes  [GiD_EntitiesGroups get $group volumes]
-    }
+    Mesh::PrepareMeshGeneration $elementsize
 
     # Change the mesh settings depending on the element requirements
-    if {[Kratos::CheckMeshCriteria $elementsize]<0} {
+    if {[Mesh::CheckMeshCriteria $elementsize]<0} {
         return "-cancel-"
     }
 
@@ -459,9 +449,17 @@ proc Kratos::Event_MeshProgress { total_percent partial_percents_0 partial_perce
 proc Kratos::Event_AfterMeshGeneration {fail} {
     variable tmp_init_mesh_time
 
+    if {$fail} {
+        Kratos::Log "Mesh generation failed"
+        return
+    }
+
     # Change the mesh settings depending on the element requirements. Reset previous settings
     # catch {Kratos::ResetMeshCriteria $fail}
 
+    Mesh::AddPointElementsIfNeeded
+
+
     # Maybe the current application needs to do some extra job
     apps::ExecuteOnCurrentApp AfterMeshGeneration $fail
     set endtime [clock seconds]

kratos.gid/scripts/Mesh.tcl

@@ -0,0 +1,111 @@
+##################################################################################
+
+namespace eval ::Mesh {
+    Kratos::AddNamespace [namespace current]
+
+
+}
+proc Mesh::PrepareMeshGeneration {elementsize} {
+
+    GiD_MeshData mesh_criteria to_be_meshed 1 lines [GiD_Geometry list line]
+    GiD_MeshData mesh_criteria to_be_meshed 1 surfaces [GiD_Geometry list surface]
+    GiD_MeshData mesh_criteria to_be_meshed 1 volumes  [GiD_Geometry list volume ]
+
+    # We need to mesh every line and surface assigned to a group that appears in the tree
+    foreach group [spdAux::GetAppliedGroups] {
+        GiD_MeshData mesh_criteria to_be_meshed 2 lines [GiD_EntitiesGroups get $group lines]
+        GiD_MeshData mesh_criteria to_be_meshed 2 surfaces [GiD_EntitiesGroups get $group surfaces]
+        GiD_MeshData mesh_criteria to_be_meshed 2 volumes  [GiD_EntitiesGroups get $group volumes]
+    }
+}
+
+
+proc Mesh::CheckMeshCriteria { elementsize } {
+
+    set force_mesh_order [dict create]
+    set elements_used [spdAux::GetUsedElements]
+    set forced_mesh_order -1
+    foreach element_id $elements_used {
+        set element [Model::getElement $element_id]
+        if {[$element hasAttribute "MeshOrder"]} {
+            set element_forces [$element getAttribute "MeshOrder"]
+            if {$element_forces eq "Quadratic"} {
+                set element_forces 1
+            } else {
+                set element_forces 0
+            }
+            dict set force_mesh_order $element_id $element_forces
+            if {$forced_mesh_order eq -1} {
+                set forced_mesh_order $element_forces
+            } else {
+                if {$forced_mesh_order ne $element_forces} {
+                    # W "The element $element_id requires a different mesh order"
+                    W "Incompatible mesh orders in elements"
+                    return -1
+                }
+            }
+        }        
+    }
+
+    if {$forced_mesh_order ne -1} {
+
+        set element [lindex [dict keys $force_mesh_order] 0]
+        set previous_mesh_order [write::isquadratic]
+        set current_mesh_type [Kratos::GetMeshOrderName $previous_mesh_order]
+        set desired_mesh_type [Kratos::GetMeshOrderName $forced_mesh_order]
+        if {$previous_mesh_order ne $forced_mesh_order} {
+            W "The element $element requires a different mesh order: $desired_mesh_type"
+            W "Currently the mesh order is $current_mesh_type. please change it in the menu Mesh > Quadratic type"
+            return -1
+        }
+    }
+    return 0
+}
+
+proc Mesh::AddPointElementsIfNeeded {} {
+    # Foreach groups assigned in tree
+    set condition_groups [spdAux::GetUsedConditions]
+
+    # condition_groups is a dict of conditionid -> list of group nodes (tdom)
+    foreach condid [dict keys $condition_groups] {
+        set cond [Model::getCondition $condid]
+        if {$cond eq ""} {
+            continue
+        }
+
+        # if group element is point and has topology for points
+        set topology [$cond getTopologyFeature "Point" 1]
+        if {$topology eq ""} {
+            continue
+        }
+
+        set group_nodes [dict get $condition_groups $condid]
+        foreach node_tdom $group_nodes {
+            set group_id [get_domnode_attribute $node_tdom n]
+
+            # if group ov is point or condition ov is point
+            set ov ""
+            if {[$node_tdom hasAttribute ov]} {
+                set ov [get_domnode_attribute $node_tdom ov]
+            } else {
+                set ov [get_domnode_attribute [$node_tdom parent] ov]
+            }
+            if { $ov ne "Point" && $ov ne "point" } {
+                continue
+            }
+
+            # Get the goodname of the group
+            set group_id [write::GetWriteGroupName $group_id]
+
+            set node_ids [GiD_EntitiesGroups get $group_id nodes]
+            set new_nodeids [list]
+            foreach nodeid $node_ids {
+                set new_nodeid [GiD_Mesh create element append Point 1 [list $nodeid]]
+                # Add to same groups as the node
+                lappend new_nodeids $new_nodeid
+            }
+            GiD_EntitiesGroups assign $group_id elements $new_nodeids
+        }
+    }
+
+}

kratos.gid/scripts/Utils.tcl

@@ -292,47 +292,6 @@ proc Kratos::GetMeshBasicData { } {
     return $result
 }
 
-proc Kratos::CheckMeshCriteria { elementsize } {
-
-    set force_mesh_order [dict create]
-    set elements_used [spdAux::GetUsedElements]
-    set forced_mesh_order -1
-    foreach element_id $elements_used {
-        set element [Model::getElement $element_id]
-        if {[$element hasAttribute "MeshOrder"]} {
-            set element_forces [$element getAttribute "MeshOrder"]
-            if {$element_forces eq "Quadratic"} {
-                set element_forces 1
-            } else {
-                set element_forces 0
-            }
-            dict set force_mesh_order $element_id $element_forces
-            if {$forced_mesh_order eq -1} {
-                set forced_mesh_order $element_forces
-            } else {
-                if {$forced_mesh_order ne $element_forces} {
-                    # W "The element $element_id requires a different mesh order"
-                    W "Incompatible mesh orders in elements"
-                    return -1
-                }
-            }
-        }        
-    }
-
-    if {$forced_mesh_order ne -1} {
-
-        set element [lindex [dict keys $force_mesh_order] 0]
-        set previous_mesh_order [write::isquadratic]
-        set current_mesh_type [Kratos::GetMeshOrderName $previous_mesh_order]
-        set desired_mesh_type [Kratos::GetMeshOrderName $forced_mesh_order]
-        if {$previous_mesh_order ne $forced_mesh_order} {
-            W "The element $element requires a different mesh order: $desired_mesh_type"
-            W "Currently the mesh order is $current_mesh_type. please change it in the menu Mesh > Quadratic type"
-            return -1
-        }
-    }
-    return 0
-}
 
 proc Kratos::GetMeshOrderName {order} {
     switch $order {

kratos.gid/scripts/Writing/WriteGeometries.tcl

@@ -40,7 +40,7 @@ proc write::printGeometryConnectivities {group etype nnodes} {
     # Prepare the indent
     set s [mdpaIndent]
     set nDim $::Model::SpatialDimension
-    set geometry_name ${etype}${nDim}${nnodes}
+    set geometry_name [GetGeometryName $etype $nDim $nnodes] 
 
     # Prepare the formats dict
     set formats [GetFormatDict $group "" $nnodes]
@@ -63,6 +63,32 @@ proc write::printGeometryConnectivities {group etype nnodes} {
         if {$etype == "Sphere" || $etype == "Circle"} {
             write::writeSphereRadiusOnGroup $group
         }
+    } else {
+        # Trick: GiD < 17.3.x return 0 if elements are of type Point
+        set elems [GiD_EntitiesGroups get $group elements -element_type point]
+        set num_elems [objarray length $elems]
+        if {$num_elems > 0} {
+            # Write header
+            WriteString "${s}Begin Geometries $geometry_name // GUI group identifier: $group"
+            # increase indent (allows folding in text editor)
+            incr ::write::current_mdpa_indent_level
+            # Write the connectivities
+            set s1 [mdpaIndent]
+            objarray foreach elem $elems {
+                set node_id [GiD_Mesh get element $elem connectivities]
+                GiD_WriteCalculationFile puts "${s1}$elem $node_id"
+            }
+            # decrease indent
+            incr ::write::current_mdpa_indent_level -1
+            # Write footer
+            WriteString "${s}End Geometries"
+            WriteString ""
+
+            # Write the radius if it is a sphere or a circle
+            if {$etype == "Sphere" || $etype == "Circle"} {
+                write::writeSphereRadiusOnGroup $group
+            }
+        }
     }
     if {[GetConfigurationAttribute time_monitor]} {set endtime [clock seconds]; set ttime [expr {$endtime-$inittime}]; W "printGeometryConnectivities $geometry_name time: [Kratos::Duration $ttime]"}
 }
@@ -73,4 +99,12 @@ proc write::writeSphereRadiusOnGroup { groupid } {
     GiD_WriteCalculationFile connectivities [dict create $groupid "%.0s %10d 0 %10g\n"]
     write::WriteString "End NodalData"
     write::WriteString ""
+}
+
+proc write::GetGeometryName { etype nDim nnodes } {
+    if {$etype == "Point"} {
+        return "${etype}${nDim}"
+    } else {
+        return "${etype}${nDim}${nnodes}"
+    }
 }

kratos.gid/scripts/spdAuxiliar.tcl

@@ -423,6 +423,29 @@ proc spdAux::GetUsedElements {{alt_un ""}} {
     return $lista
 }
 
+# returns ["condition1" {gNode1 gNode2} "condition2" {gNode3 gNode4}]
+proc spdAux::GetUsedConditions {{root ""}} {
+    set resultDict [dict create]
+
+
+    set xp1 "./condition/group"
+    if {$root eq "" } {
+        set root [customlib::GetBaseRoot]
+        set xp1 "//condition/group"
+    }
+
+    foreach gNode [$root selectNodes $xp1] {
+        set condition_node [$gNode parent]
+        set condition_name [$condition_node @n]
+
+        set cond [Model::getCondition $condition_name]
+        if {$cond eq ""} {continue}
+        dict lappend resultDict $condition_name $gNode
+
+    }
+    return $resultDict
+}
+
 proc spdAux::LoadIntervalGroups { {root ""} } {
     customlib::UpdateDocument
     variable GroupsEdited