Support for Kratos mdpa files

pygem/__init__.py

@@ -20,4 +20,5 @@
 from .stephandler import StepHandler
 from .igeshandler import IgesHandler
 from .khandler import KHandler
+from .mdpahandler import MdpaHandler
 from .params import *

pygem/mdpahandler.py

@@ -0,0 +1,90 @@
+"""
+Derived module from khandler.py to handle Kratos mesh (.mdpa) files.
+"""
+import numpy as np
+import pygem.filehandler as fh
+
+
+class MdpaHandler(fh.FileHandler):
+    """
+    Kratos mesh file handler class
+
+    :cvar string infile: name of the input file to be processed.
+    :cvar string outfile: name of the output file where to write in.
+    :cvar list extensions: extensions of the input/output files. It is equal
+            to '.mdpa'.
+    """
+
+    def __init__(self):
+        super(MdpaHandler, self).__init__()
+        self.extensions = ['.mdpa']
+
+    def parse(self, filename):
+        """
+        Method to parse the file `filename`. It returns a matrix with all the
+        coordinates. It reads only the section after "Begin Nodes" of the
+        mdpa files.
+
+        :param string filename: name of the input file.
+
+        :return: mesh_points: it is a `n_points`-by-3 matrix containing the
+                coordinates of the points of the mesh.
+        :rtype: numpy.ndarray
+        """
+        self._check_filename_type(filename)
+        self._check_extension(filename)
+        self.infile = filename
+        index = -9
+        mesh_points = []
+        with open(self.infile, 'r') as input_file:
+            for num, line in enumerate(input_file):
+                if line.startswith('Begin Nodes'):
+                    index = num
+                if num == index + 1:
+                    if line.startswith('End Nodes'):
+                        break
+                    else:
+                        line = line.replace('D', 'E')
+                        li = []
+                        for t in line.split()[1:]:
+                            try:
+                                li.append(float(t))
+                            except ValueError:
+                                pass
+                        mesh_points.append(li)
+                        index = num
+            mesh_points = np.array(mesh_points)
+        return mesh_points
+
+    def write(self, mesh_points, filename):
+        """
+        Writes a .mdpa file, called filename, copying all the lines from
+        self.filename but the coordinates. mesh_points is a matrix that
+        contains the new coordinates to write in the .mdpa file.
+
+        :param numpy.ndarray mesh_points: it is a `n_points`-by-3 matrix
+                          containing the coordinates of the points of the mesh
+        :param string filename: name of the output file.
+        """
+        self._check_filename_type(filename)
+        self._check_extension(filename)
+        self._check_infile_instantiation()
+        self.outfile = filename
+        index = -9
+        i = 0
+        with open(self.outfile, 'w') as output_file:
+            with open(self.infile, 'r') as input_file:
+                for num, line in enumerate(input_file):
+                    if line.startswith('Begin Nodes'):
+                        index = num
+                    if num == index + 1:
+                        if line.startswith('End Nodes'):
+                            index = -9
+                        else:
+                            line = (" {:6d} {:23.16E} {:23.16E} {:23.16E}\n"
+                                    .format(i+1, mesh_points[i][0],
+                                            mesh_points[i][1],
+                                            mesh_points[i][2]))
+                            i += 1
+                            index = num
+                    output_file.write(line)

tests/test_mdpahandler.py

@@ -0,0 +1,122 @@
+from unittest import TestCase
+import pygem.mdpahandler as uh
+import numpy as np
+import filecmp
+import os
+
+
+class TestMdpaHandler(TestCase):
+    def test_mdpa_instantiation(self):
+        mdpa_handler = uh.MdpaHandler()
+
+    def test_mdpa_default_infile_member(self):
+        mdpa_handler = uh.MdpaHandler()
+        self.assertIsNone(mdpa_handler.infile)
+
+    def test_mdpa_default_outfile_member(self):
+        mdpa_handler = uh.MdpaHandler()
+        self.assertIsNone(mdpa_handler.outfile)
+
+    def test_mdpa_default_extension_member(self):
+        mdpa_handler = uh.MdpaHandler()
+        self.assertListEqual(mdpa_handler.extensions, ['.mdpa'])
+
+    def test_mdpa_parse_failing_filename_type(self):
+        mdpa_handler = uh.MdpaHandler()
+        with self.assertRaises(TypeError):
+            mesh_points = mdpa_handler.parse(5.2)
+
+    def test_mdpa_parse_failing_check_extension(self):
+        mdpa_handler = uh.MdpaHandler()
+        with self.assertRaises(ValueError):
+            mesh_points = mdpa_handler.parse(
+                'tests/test_datasets/test_square.iges')
+
+    def test_mdpa_parse_infile(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        self.assertEqual(mdpa_handler.infile, 'tests/test_datasets/test_square.mdpa')
+
+    def test_mdpa_parse_shape(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        self.assertTupleEqual(mesh_points.shape, (256, 3))
+
+    def test_mdpa_parse_coords_1(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        np.testing.assert_almost_equal(mesh_points[190][0], 1.0)
+
+    def test_mdpa_parse_coords_2(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        np.testing.assert_almost_equal(mesh_points[72][1], 0.4)
+
+    def test_mdpa_parse_coords_3(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        np.testing.assert_almost_equal(mesh_points[100][2], 0.0)
+
+    def test_mdpa_parse_coords_4(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        np.testing.assert_almost_equal(mesh_points[0][0], 0.0)
+
+    def test_mdpa_parse_coords_5(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        np.testing.assert_almost_equal(mesh_points[-1][2], 0.0)
+
+    def test_mdpa_write_failing_filename_type(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        with self.assertRaises(TypeError):
+            mdpa_handler.write(mesh_points, -2)
+
+    def test_mdpa_write_failing_check_extension(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse('tests/test_datasets/test_square.mdpa')
+        with self.assertRaises(ValueError):
+            mdpa_handler.write(mesh_points, 'tests/test_datasets/test_square.iges')
+
+    def test_mdpa_write_failing_infile_instantiation(self):
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = np.zeros((20, 3))
+        with self.assertRaises(RuntimeError):
+            mdpa_handler.write(mesh_points,
+                            'tests/test_datasets/test_square_out.mdpa')
+
+    def test_mdpa_write_outfile(self):
+        infilename = 'tests/test_datasets/test_square.mdpa'
+        outfilename = 'tests/test_datasets/test_square_out.mdpa'
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse(infilename)
+        mdpa_handler.write(mesh_points, outfilename)
+        self.assertEqual(mdpa_handler.outfile, outfilename)
+        #self.addCleanup(os.remove, outfilename)
+
+    def test_mdpa_write_comparison_1(self):
+        infilename = 'tests/test_datasets/test_square.mdpa'
+        outfilename = 'tests/test_datasets/test_square_out.mdpa'
+        outfilename_expected = 'tests/test_datasets/test_square.mdpa'
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse(infilename)
+        mdpa_handler.write(mesh_points, outfilename)
+        self.assertTrue(filecmp.cmp(outfilename, outfilename_expected))
+        self.addCleanup(os.remove, outfilename)
+
+    def test_mdpa_write_comparison_2(self):
+        infilename = 'tests/test_datasets/test_square.mdpa'
+        outfilename = 'tests/test_datasets/test_square_out.mdpa'
+        outfilename_expected = 'tests/test_datasets/test_square_out_true.mdpa'
+        mdpa_handler = uh.MdpaHandler()
+        mesh_points = mdpa_handler.parse(infilename)
+        mesh_points[0][0] = 0.0
+        mesh_points[5][1] = 1.0
+        mesh_points[9][2] = 0.0
+        mesh_points[44][0] = 0.0
+        mesh_points[122][1] = 0.2
+        mesh_points[255][2] = 0.0
+        mdpa_handler.write(mesh_points, outfilename)
+        self.assertTrue(filecmp.cmp(outfilename, outfilename_expected))
+        self.addCleanup(os.remove, outfilename)

tutorials/README.md

@@ -17,6 +17,11 @@ Here we present the radial basis functions interpolation technique on a simple c
 #### [Tutorial 5](https://github.com/mathLab/PyGeM/blob/master/tutorials/tutorial-5-idw.ipynb)
 Here we show the inverse distance weighting interpolation technique on a simple cube.
 
+#### [Tutorial 6](https://github.com/mathLab/PyGeM/blob/master/tutorials/tutorial-6-k.ipynb)
+This tutorial shows how to deform a LS-Dyna k file with a prescribed continuity using the free-form deformation.
+
+#### [Tutorial 7](https://github.com/mathLab/PyGeM/blob/master/tutorials/tutorial-7-mdpa.ipynb)
+This tutorial shows how to deform a Kratos Multiphysics mdpa file with a prescribed continuity using the free-form deformation. It is shown how these can then be visualized using the Python library VTK Plotter.
 
 #### More to come...
 We plan to add more tutorials but the time is often against us. If you want to contribute with a notebook on a feature not covered yet we will be very happy and give you support on editing!