Tutorial: Finite Element Creation¶

The Finite Element Creator application is designed to illustrate how to create a 2D finite element and visualise it.

This tutorial shows how to

create a 2D finite element using a template
visualise the element using a surface graphic and node points

The souce code used in this tutorial is available from the physiome project svn server.

../../_images/finite_element_creation.png

Overview¶

Finite elements are the building blocks of any model. To create a model from finite elements the Zinc library uses templates. Using templates enables a modeller to create very large model with ease.

This tutorial uses the same technique for creating an application as the simple axis viewer tutorial, see this tutorial for more information.

Initialise¶

The initialisation is done in two parts one for the PyZinc context and one for the OpenGL Context. Creating the PyZinc context is straightforward, we need only create a Context object and give it a unique name:

self._context = Context("finiteelementcreation")

We keep a handle to the PyZinc context as any handles originating from this context will be invalid once the context has been deleted.

The OpenGL Context is created by Qt and we must set any Zinc scene viewers that we create to have the same properties.

    def initializeGL(self):
        '''
        Initialise the Zinc scene, create the finite element, and the surface to visualise it.  
        '''
               
        # From the graphics module get the scene viewer module.
        scene_viewer_module = self._context.getSceneviewermodule()
        
        # Get the glyph module from the graphics module and define the standard glyphs
        glyph_module = self._context.getGlyphmodule()
        glyph_module.defineStandardGlyphs()

        # From the scene viewer package we can create a scene viewer, we set up the scene viewer to have the same OpenGL properties as
        # the QGLWidget.
        self._scene_viewer = scene_viewer_module.createSceneviewer(Sceneviewer.BUFFERING_MODE_DOUBLE, Sceneviewer.STEREO_MODE_DEFAULT)

        # Create a filter for visibility flags which will allow us to see our graphics.
        filter_module = self._context.getScenefiltermodule()
        graphics_filter = filter_module.createScenefilterVisibilityFlags()

        # Set the filter for the scene viewer to use.
        self._scene_viewer.setScenefilter(graphics_filter)
        
        self.create2DFiniteElement()
        self.createSurfaceGraphic()

        # Place the viewport to contain everything visible in the scene.
        self._scene_viewer.viewAll()
        self._notifier = self._scene_viewer.createSceneviewernotifier()
        self._notifier.setCallback(self.sceneviewerEvent)

Line 10 we get the scene viewer package from the graphics module, from the scene viewer package we can create scene viewers.

Line 18 here we create a scene viewer we also keep the handle that is returned in the class because we will use it a lot.

Line 25 function call to create a basic linear 2D finite element

Line 26 function call to create a surface graphic so that we may visualise the 2D finite element created on line 25.

Line 29 using the viewAll() API function from the scene viewer we can set the viewport to contain all the visible elements of the scene viewer’s scene.

Create¶

To create a finite element we must use templates, this is a bit cumbersome for creating basic meshes but it is much more useful when we create very big meshes. Here is the code:

    def create2DFiniteElement(self):
        # Get a the root region, which is the default region of the context.
        default_region = self._context.getDefaultRegion()
        
        # Get the field module for root region, with which we  shall create a 
        # finite element coordinate field.
        field_module = default_region.getFieldmodule()
        
        field_module.beginChange()
        
        # Create a finite element field with 2 components to represent 2 dimensions
        finite_element_field = field_module.createFieldFiniteElement(2)
        # Set the name of the field, we give it label to help us understand it's purpose
        finite_element_field.setName('coordinates')
        self._finite_element_field = finite_element_field
        # Find a special node set named 'cmiss_nodes'
        nodeset = field_module.findNodesetByName('nodes')
        node_template = nodeset.createNodetemplate()
        # Set the finite element coordinate field for the nodes to use
        node_template.defineField(finite_element_field)
        
        field_cache = field_module.createFieldcache()
        
        # The node template for the square expects the nodes to be ordered in 
        # a mirrored Z - pattern for the square element.  The order of the coordinates
        # for the nodes is unimportant here but it will allow us to index the nodes in
        # order later when using the element template.
        node_coordinates = [[0, 0.0], [3.0, 0.0], [0.0, 4.0], [2.0, 2.0]]
        # Create four nodes to define a square finite element over
        for i, node_coordinate in enumerate(node_coordinates):
            # Node indexes start from 1
            node = nodeset.createNode(i+1, node_template)
            # Set the node coordinates, first set the field cache to use the current node
            field_cache.setNode(node)
            # Pass in floats as an array
            finite_element_field.assignReal(field_cache, node_coordinate)

        # We want to create a 2D element so we grab the predefined 2D mesh.  Currently there
        # is only one mesh for each dimension 1D, 2D, and 3D the user is not able to name
        # their own mesh.  This is due to change in an upcoming release of PyZinc.
        mesh = field_module.findMeshByDimension(2)
        element_template = mesh.createElementtemplate()
        element_template.setElementShapeType(Element.SHAPE_TYPE_SQUARE)
        element_node_count = 4
        element_template.setNumberOfNodes(element_node_count)
        # Specify the dimension and the interpolation function for the element basis function. 
        linear_basis = field_module.createElementbasis(2, Elementbasis.FUNCTION_TYPE_LINEAR_LAGRANGE)
        # The indexes of the nodes in the node template we want to use
        node_indexes = [1, 2, 3, 4]
        # Define a nodally interpolated element field or field component in the
        # element_template. Only Lagrange, simplex and constant basis function types
        # may be used with this function, i.e. where only a simple node value is
        # mapped. Shape must be set before calling this function.  The -1 for the component number
        # defines all components with identical basis and nodal mappings.
        element_template.defineFieldSimpleNodal(finite_element_field, -1, linear_basis, node_indexes)
                    
        for i in range(element_node_count):
            node = nodeset.findNodeByIdentifier(i+1)
            element_template.setNode(i+1, node)

        mesh.defineElement(-1, element_template)
        
        field_module.endChange()

We will cover some of the more pertinent parts of this code rather taking it line by line.

Line 12 from the field module we create a finite element field with the desired number of components. For us the number of components will be used to give us values for a two-dimensional coordinate system.

Line 14 set the name of the finite element field to ‘coordinates’. We can then search the field module for this name later on when we require it.

Line 16 by setting the is managed flag we are telling the field module to keep this field around for us and not delete it when it holds the only reference.

Lines 19, 20 there exists a special nodeset named ‘cmiss_nodes’ where we can create nodes that can used for constructing finite elements. From the nodeset we create a node template for creating the nodes required for our finite element.

Line 23 we set the finite element field as the field for the nodes to use. We can add any number of fields to the node template at this point.

Line 25 create a field cache for storing or retrieving known field values at a location.

Line 31 define the coordinates of the nodes for our finite element. Because we are goint to define a square two-dimensional finite element with linear basis functions we will use four distinct nodes. We have a range of options open to us here repeated nodes, quadratic basis functions etc. but in this situation we will use four distinct nodes to describe our finite element.

Lines 33 - 39 we create our nodes from the nodeset by giving them a unique identifier in the set and the template to create them from. We also set the value of the field at the node location by first setting the node into the field cache and then assigning the value of the node coordinate to all items in the field cache. In this case this is the current node we have just created.

Lines 44 - 48 we want to create a top-level 2D element so we get the two-dimensional mesh. We could create a top-level 3D or 1D element similarly. From the mesh we create an element template and set the element shape to square. We also tell the template how many nodes we are going to use for defining elements created from the template.

Line 50 create the element basis function which we are going to use for elements created from this template.

Line 52 define the indexes of the nodes in the element template in, the node index order is important to create an element. The nodes must be given so that the lowest xi coordinate varies fastest. For the square with linear basis functions we have bottom-left, bottom-right, top-left, top-right ordering.

Line 58 define a simple nodal field on the element template over the finite element field for all components using a linear basis with node indexes [1, 2, 3, 4].

Lines 60 - 63 we set the nodes from the nodeset into the element template at the given node index.

Line 64 define the element automatically generating the integer identifier from the element template. If we required a handle to the element we could use the createElement() function from the mesh object.

Visualise¶

Once we have created a finite element we would like to visualise it. Because we have made a 2D finite element we can view it using a surface graphic. Much of this code is the same as we have seen in the simple axis viewer tutorial so we will only look at the bits that have changed.

        field_module = default_region.getFieldmodule()
        finite_element_field = field_module.findFieldByName('coordinates')
        # Create a surface graphic and set it's coordinate field to the finite element coordinate field
        # named coordinates
        surface = scene.createGraphicsSurfaces()
        surface.setCoordinateField(finite_element_field)
        
        # Create point graphics and set the coordinate field to the finite element coordinate field
        # named coordinates
        sphere = scene.createGraphicsPoints()
        sphere.setCoordinateField(finite_element_field)
        sphere.setFieldDomainType(Field.DOMAIN_TYPE_NODES)
        att = sphere.getGraphicspointattributes()
        att.setGlyphShapeType(Glyph.SHAPE_TYPE_SPHERE)
        att.setBaseSize([1])

Line 1 get the field module for the region.

Line 2 we get the finite element field created in this field module by name. Alternatively we could have kept a handle to the field created in create2DFiniteElement().

Line 6 from the scene for the region create a surface graphic.

Line 7 set the coordinate field for the graphic. The finite element field is defined over the domain of the finite element, the domain is specified by the nodes we created.

Line 11 from the scene create a graphic points visualisation.

Line 13 set the domain type for the graphic points to be DOMAIN_NODES.

Line 15 set the glyph type to GLYPH_TYPE_SPHERE, note: if we don’t call defineStandardGlyphs() then we won’t see the node points.

Line 16 set the base size of the graphic, again note: if we don’t set the base size we won’t see the node points