0404

msgpi/ms/beam.py

@@ -1,6 +1,8 @@
 import numpy as np
 import msgpi.sg as sgi
-import msgpi.io.utils as utl
+import msgpi.io.utils as ioutil
+import msgpi.utils as utils
+import pprint
 
 
 class BeamSegment(object):
@@ -18,6 +20,8 @@ def __init__(self):
         #: `[starting, ending]`
         self.css = []
 
+        self.cs_template = []
+
         #: float: Rotation around a1.
         self.rotate_a1 = 0.0
         #: float: Twist.
@@ -34,6 +38,15 @@ def __init__(self):
         #: int: Number of division of the segment.
         self.num_divisions = 1
 
+        self.cs_parameters = {}  #: Parameter distributions
+
+    def summary(self):
+        print('  points:', self.points)
+        print('  coords:', self.coords)
+        print('  cs names:', self.css)
+        print('  cs templates:', self.cs_template)
+        print('  cs parameters:', self.cs_parameters)
+
     def calcLengthSq(self):
         """Calculate the square of the segment length.
 
@@ -96,9 +109,34 @@ def __init__(self):
         self.timefunctions = []  #: Time functions
         self.initcurvatures = []  #: Initial curvatures
 
+        self.functions = {}  #: Functions
+        self.distributions = {}  #: Distributions
+
         # Results
         # self.results = None  #: Results of GEBT analysis
 
+
+    def summary(self):
+        pp = pprint.PrettyPrinter(indent=1, width=80)
+
+        print('\nPoints')
+        pp.pprint(self.points)
+
+        print('\nDistributions')
+        pp.pprint(self.distributions)
+
+        print('\nFunctions')
+        # pp.pprint(self.functions)
+        for k, v in self.functions.items():
+            print(k, ':', v)
+
+        print('\nSegments')
+        for k, v in self.segments.items():
+            print('segment', k)
+            v.summary()
+
+
+
     def echo(self):
         """Print the beam data.
         """
@@ -179,15 +217,15 @@ def writeGmshMsh(self):
             # Reference line
             for i, st in enumerate(self.stations):
                 fout.write('{0:8d}'.format(i+1))
-                utl.writeFormatFloats(fout, np.array(st['coordinates']))
+                ioutil.writeFormatFloats(fout, np.array(st['coordinates']))
             # Cross sections
             nnodes = len(self.stations)
             for st in self.stations:
                 loc = np.array(st['coordinates'])
                 cs = st['cross_section']
                 for nid, coords in cs.nodes.items():
                     fout.write('{0:8d}'.format(nid+nnodes))
-                    utl.writeFormatFloats(fout, np.array(coords)+loc)
+                    ioutil.writeFormatFloats(fout, np.array(coords)+loc)
                 nnodes += cs.num_nodes
             fout.write('$EndNodes\n')
             fout.write('\n')
@@ -198,7 +236,7 @@ def writeGmshMsh(self):
             # Reference line
             for i in range(len(self.stations) - 1):
                 line = [i+1, 1, 2, 1, 1, i+1, i+2]
-                utl.writeFormatIntegers(fout, line)
+                ioutil.writeFormatIntegers(fout, line)
             # Cross sections
             nnodes = len(self.stations)
             nelems = len(self.stations) - 1
@@ -207,7 +245,7 @@ def writeGmshMsh(self):
                 for eid, nodes in cs.elements.items():
                     line = np.array([eid+nelems, 2, 2, 1, 1])
                     line = np.concatenate([line, np.array(nodes)+nnodes])
-                    utl.writeFormatIntegers(fout, line)
+                    ioutil.writeFormatIntegers(fout, line)
                 nnodes += cs.num_nodes
                 nelems += cs.num_elements
             fout.write('$EndElements\n')

msgpi/ms/prebeam.py

@@ -4,7 +4,7 @@
 import msgpi.ms.beam as msb
 import msgpi.sg as sgi
 import msgpi.cross_section as sgcs
-import msgpi.utils as utl
+import msgpi.utils as utils
 
 
 def preBeam(fn_beam, mode=1, sections=[]):
@@ -477,3 +477,130 @@ def preBeam(fn_beam, mode=1, sections=[]):
     #     self.divisions.append(int(np.ceil(dx1 / ms)))
 
     return beam
+
+
+
+
+
+
+
+
+
+def preBeam2(fn_beam):
+
+    beam = msb.Beam()
+
+    with open(fn_beam, 'r') as fo:
+        tree = et.parse(fo)
+        xr_beam = tree.getroot()
+
+        # Read design
+        xe_design = xr_beam.find('design')
+
+
+
+
+        # Read points
+        # pn2l = {}  # point name to label dictionary
+        # plabel = 0
+        xe_points = xe_design.find('points')
+        for xe_point in xe_points.findall('point'):
+            # plabel += 1
+            pname = xe_point.get('name')
+            pcoord = list(map(float, xe_point.text.strip().split()))
+            # pn2l[pname] = plabel
+            # beam.points[plabel] = pcoord
+            beam.points[pname] = pcoord
+
+
+
+
+        # Read functions
+        xe_functions = xr_beam.find('functions')
+        for xe_function in xe_functions.findall('function'):
+            fname = xe_function.get('name')
+            ftype = xe_function.get('type', default='polynomial')
+            fpointsx = []
+            fpointsy = []
+
+            xe_coeffs = xe_function.find('coefficients')
+            if xe_coeffs is not None:
+                # print(xe_coeffs.text)
+                fcoeffs = list(map(float, xe_coeffs.text.strip().split()))
+            xe_points = xe_function.find('points')
+            if xe_points is not None:
+                for xe_point in xe_points.findall('point'):
+                    x, y = list(map(float, xe_point.text.strip().split()))
+                    fpointsx.append(x)
+                    fpointsy.append(y)
+
+            beam.functions[fname] = utils.PolynomialFunction(
+                coefficients=fcoeffs, x=fpointsx, y=fpointsy
+            )
+
+
+
+
+        # Read distributions
+        xe_distributions = xe_design.find('distributions')
+        for xe_distribution in xe_distributions.findall('distribution'):
+            # distr_type = xe_distribution.get('type', default='function')
+            param_name = xe_distribution.find('parameter').text
+
+            func_name = xe_distribution.find('function').text
+            m = 1.0
+            if xe_distribution.find('multiplier'):
+                m = float(xe_distribution.find('multiplier').text)
+            beam.distributions[param_name] = [func_name, m]
+
+
+
+
+        # Read segments
+        xe_segments = xe_design.find('segments')
+        bslabel = 0
+        for xe_sgm in xe_segments.findall('segment'):
+            bslabel += 1
+            bs = msb.BeamSegment()
+
+            # Beginning station
+            xe_start = xe_sgm.find('begin')
+            pn = xe_start.find('location').text.strip()
+            # bs.points.append(pn2l[pn])
+            bs.points.append(pn)
+            # pc = beam.findPtCoordByName(pn2l[pn])
+            pc = beam.findPtCoordByName(pn)
+            bs.coords.append(pc)
+            xe_cs = xe_start.find('cross_section')
+            fn_cs_template = xe_cs.get('template')
+            cs_name = xe_cs.text.strip()
+            # fn_cs = cs_name + '.xml'
+            bs.cs_template.append(fn_cs_template)
+            bs.css.append(cs_name)
+
+            # Ending station
+            xe_stop = xe_sgm.find('end')
+            pn = xe_stop.find('location').text.strip()
+            # bs.points.append(pn2l[pn])
+            bs.points.append(pn)
+            # pc = beam.findPtCoordByName(pn2l[pn])
+            pc = beam.findPtCoordByName(pn)
+            bs.coords.append(pc)
+            xe_cs = xe_stop.find('cross_section')
+            fn_cs_template = xe_cs.get('template')
+            cs_name = xe_cs.text.strip()
+            # fn_cs = cs_name + '.xml'
+            bs.cs_template.append(fn_cs_template)
+            bs.css.append(cs_name)
+
+            # Cross-section parameters
+            for xe_param in xe_sgm.findall('parameter'):
+                name = xe_param.find('name').text.strip()
+                value = float(xe_param.find('value').text.strip())
+                bs.cs_parameters[name] = value
+
+
+            beam.segments[bslabel] = bs
+
+
+    return beam

msgpi/utils.py

@@ -336,3 +336,43 @@ def __call__(self, x):
         basis = np.asarray(basis)
 
         return np.dot(self.coeffs, basis)
+
+
+
+
+
+
+
+
+
+class PolynomialFunction(object):
+    def __init__(self, coefficients=[], x=[], y=[]):
+        self.coefficients = coefficients
+        self.points = {}
+        if len(x) > 0:
+            # points = np.asarray(points)
+            ndof = len(x)
+            A = np.zeros((ndof, ndof))
+            for i in range(ndof):
+                self.points[x[i]] = y[i]
+                for j in range(ndof):
+                    A[i, j] = x[i]**j
+            # print(A)
+            self.coefficients = np.linalg.solve(A, y)
+            # print(self.coefficients)
+
+    def __call__(self, x):
+        f = 0.0
+        if x in self.points.keys():
+            f = self.points[x]
+        else:
+            for i, c in enumerate(self.coefficients):
+                f += c * x**i
+        return f
+
+    def __str__(self):
+        terms = []
+        for i, c in enumerate(self.coefficients):
+            terms.append(f'({c})*x^{i}')
+
+        return ' + '.join(terms)