Openmc alpha

GDML_math.md

@@ -0,0 +1,61 @@
+In GDML (Geometry Description Markup Language), which is used for defining complex geometries in simulations, particularly within the context of particle physics (like with Geant4), certain mathematical functions are allowed for specifying parameters and dimensions. These functions are used in expressions within GDML to provide dynamic values for various geometry properties.
+
+### Allowed Mathematical Functions in GDML:
+GDML supports a limited set of standard mathematical functions that you can use directly within your XML definitions. These functions include:
+
+1. **Arithmetic Operations:**
+   - `+` : Addition
+   - `-` : Subtraction
+   - `*` : Multiplication
+   - `/` : Division
+   - `%` : Modulus (remainder)
+
+2. **Power Functions:**
+   - `pow(a, b)` : Returns `a` raised to the power of `b`.
+
+3. **Square Root:**
+   - `sqrt(x)` : Returns the square root of `x`.
+
+4. **Trigonometric Functions:**
+   - `sin(x)` : Sine of `x` (angle in radians).
+   - `cos(x)` : Cosine of `x` (angle in radians).
+   - `tan(x)` : Tangent of `x` (angle in radians).
+   - `asin(x)` : Arcsine of `x`, result in radians.
+   - `acos(x)` : Arccosine of `x`, result in radians.
+   - `atan(x)` : Arctangent of `x`, result in radians.
+   - `atan2(y, x)` : Arctangent of `y/x`, result in radians (takes into account the signs of both arguments to determine the correct quadrant).
+
+5. **Hyperbolic Functions:**
+   - `sinh(x)` : Hyperbolic sine of `x`.
+   - `cosh(x)` : Hyperbolic cosine of `x`.
+   - `tanh(x)` : Hyperbolic tangent of `x`.
+
+6. **Exponential and Logarithmic Functions:**
+   - `exp(x)` : Returns `e` raised to the power of `x`.
+   - `log(x)` : Natural logarithm of `x` (log base `e`).
+   - `log10(x)` : Base-10 logarithm of `x`.
+
+7. **Absolute Value:**
+   - `abs(x)` : Returns the absolute value of `x`.
+
+8. **Minimum and Maximum:**
+   - `min(a, b)` : Returns the smaller of `a` and `b`.
+   - `max(a, b)` : Returns the larger of `a` and `b`.
+
+9. **Constants:**
+   - `pi` : Represents the mathematical constant π (pi).
+   - `e` : Represents the base of the natural logarithm (approximately 2.71828).
+
+### Usage in GDML:
+These functions can be used in attributes for defining dimensions, positions, rotations, and other numeric properties. For example:
+
+```xml
+<box name="MyBox" x="2*sqrt(2)" y="5" z="10"/>
+<rotation name="MyRotation" x="cos(pi/4)" y="sin(pi/4)" z="0"/>
+```
+
+### Notes:
+- **Syntax:** Functions should be written in a syntax that matches typical programming or mathematical conventions (e.g., `sqrt(4)` for the square root of 4).
+- **Units:** When using these functions, ensure that any numerical results are in the correct units, as GDML does not automatically handle unit conversions.
+
+These functions give GDML users flexibility in defining complex geometrical shapes and relationships, ensuring that geometrical parameters can be dynamically calculated based on expressions rather than fixed values.

Utils/calcCenterOfMass.py

@@ -166,7 +166,7 @@ def MonteCarloMoments(shape, nsim: int) -> Moments:
     return moments
 
 
-def topObjects(part: App.Part) -> list:
+def topObjects(part) -> list:
     ''' return a list containing the top shapes contained in part
     a top shape may contain other shapes, but cannot be contained in
     other shapes within Part
@@ -177,7 +177,7 @@ def topObjects(part: App.Part) -> list:
         return []
 
 
-def partCM(part: App.Part) -> tuple:
+def partCM(part) -> tuple:
     '''
     Note, what we are calling the total moment of inertia is NOT the total moment
     of inertia at all!
@@ -251,6 +251,26 @@ def partCM(part: App.Part) -> tuple:
 
             # prettyPrint(obj.Label, vol, cmglob, IIorigin)
 
+        elif obj.TypeId == "Mesh::Feature":
+            mesh = obj.Mesh
+            shape = Part.Shape()
+            shape.makeShapeFromMesh(mesh.Topology, 0.5)
+            solid = Part.makeSolid(shape)
+            II0 = solid.MatrixOfInertia
+            globalPlacement = obj.getGlobalPlacement()
+            rotglob = globalPlacement.Rotation
+            rotobj = obj.Placement.Rotation
+            cm0 = shape.CenterOfGravity
+            vol = shape.Volume
+            cmglob = globalPlacement.Base + rotglob*rotobj.inverted()*(cm0 - obj.Placement.Base)
+            mom0: Moments = Moments(vol, cmglob, II0)
+            mom0.M = mom0.inertiaRotated(rotobj.inverted())
+            mom0.M = mom0.inertiaRotated(rotglob)
+            IIorigin: Matrix = mom0.inertiaAboutOrigin()
+            II += IIorigin
+            cm += vol * cmglob
+            totVol += vol
+
         elif obj.TypeId == 'App::Part':
             partVol, cm0, II0 = partCM(obj)
             cm += partVol * cm0

freecad/gdml/GDMLMaterials.py

@@ -31,6 +31,8 @@
 
 from PySide import QtGui, QtCore
 
+from freecad.gdml.importGDML import processReactor
+
 
 class GDMLMaterial(QtGui.QComboBox):
 
@@ -62,10 +64,22 @@ def getMaterialsList():
         materials = doc.Materials
         g4Mats = doc.getObject('G4Materials')
 
+    try:
+        reactorMats = doc.getObject('ReactorMaterials')
+
+    except:
+        from .importGDML import processReactor
+        from .init_gui import joinDir
+
+        print('Load Geant4 Materials XML')
+        processReactor(doc, joinDir("Resources/ReactorMaterials.xml"))
+        materials = doc.Materials
+        reactorMats = doc.getObject('ReactorMaterials')
+
     try:
         if materials is not None:
             for m in materials.OutList:
-                if m.Label != "Geant4":
+                if m.Label != "Geant4" and m.Label != "RaactorMaterials":
                     matList.append(m.Label)
                     print(matList)
     except:
@@ -80,6 +94,15 @@ def getMaterialsList():
     except:
         pass
 
+    try:
+        if reactorMats is not None:
+            for m in reactorMats.OutList:
+                for n in m.OutList:
+                    matList.append(n.Label)
+            # print(matList)
+    except:
+        pass
+
     return matList
 
 
@@ -100,20 +123,23 @@ def refreshG4Materials(doc):
 
 
 def newGetGroupedMaterials():
-    from .importGDML import joinDir, processGEANT4
+    from .importGDML import joinDir, processGEANT4, processReactor
     from .GDMLObjects import GroupedMaterials
     print(f'New getGroupedMaterials len GroupMaterials {len(GroupedMaterials)}')
     # if len(GroupedMaterials) == 0:
     mlen = len(GroupedMaterials)
     if mlen >= 0:
         doc = FreeCAD.activeDocument()
-        if not hasattr(doc, 'Materials') or not hasattr(doc, 'G4Materials'):
+        if not hasattr(doc, 'Materials') or not hasattr(doc, 'G4Materials') or not hasattr(doc, 'ReactorMaterials'):
             processGEANT4(doc, joinDir("Resources/Geant4Materials.xml"))
+            processReactor(doc, joinDir("Resources/ReactorMaterials.xml"))
             docG4Materials = doc.G4Materials
             if not hasattr(docG4Materials, 'version'):
                 refreshG4Materials(doc)
         docG4Materials = doc.G4Materials
+        reactorMaterials = doc.ReactorMaterials
         print(f'doc.G4Materials {docG4Materials}')
+        print(f'doc.ReactorMaterials {reactorMaterials}')
         for g in docG4Materials.Group:
             # print(f'g : {g.Label}')
             for s in g.Group:
@@ -122,6 +148,9 @@ def newGetGroupedMaterials():
                     GroupedMaterials[g.Label].append(s.Label)
                 else:
                     GroupedMaterials[g.Label] = [s.Label]
+
+        GroupedMaterials[reactorMaterials.Label] = [g.Label for g in reactorMaterials.Group]
+
         matList = []
         docMaterials = doc.Materials
         print(f'doc.Materials {docMaterials}')
@@ -136,44 +165,3 @@ def newGetGroupedMaterials():
             GroupedMaterials['Normal'] = matList
 
     return GroupedMaterials
-
-
-def getGroupedMaterials():
-    print('getGroupedMaterials')
-    from .GDMLObjects import GroupedMaterials
-    from .importGDML import setupEtree
-    from .init_gui import joinDir
-
-    if len(GroupedMaterials) == 0:
-        etree, root = setupEtree(joinDir("Resources/Geant4Materials.xml"))
-        materials = root.find('materials')
-
-        for material in materials.findall('material'):
-            name = material.get('name')
-            print(name)
-            if name is None:
-                print("Missing Name")
-            else:
-                for auxiliary in material.findall('auxiliary'):
-                    auxtype = auxiliary.get('auxtype')
-                    if auxtype == 'Material-type':
-                        auxvalue = auxiliary.get('auxvalue')
-                        if auxvalue in GroupedMaterials:
-                            GroupedMaterials[auxvalue].append(name)
-                        else:
-                            GroupedMaterials[auxvalue] = [name]
-
-    doc = FreeCAD.activeDocument()
-    docMaterials = doc.Materials
-    matList = []
-
-    if doc.Materials is not None:
-        for m in docMaterials.OutList:
-            if m.Label != "Geant4":
-                if m.Label not in matList:
-                    matList.append(m.Label)
-
-    if len(matList) > 0:
-        GroupedMaterials['Normal'] = matList
-
-    return GroupedMaterials