Opencmiss mesh solver for a defined direction

moveNodeToMesh.py

    def moveNode(region, nodekey, plane_normal=[0, 1, 0], cache=none):
        # moveNode uses dot products combined with opencmiss' evaluateMeshLocation function to solve where a point lies along a given normal (line in 3D).
      
        # usage: Please not that this solver assumes that the user knows where the solution should be in one dimension, 
        # for example: Where is the closest mesh point at x=0 (normal=[1,0,0]) starting at [1,3,2]? 
        
        # Inputs: 
        #   region: the region you wish to solve in (must contain a mesh and a the node you wish to move)
        #   nodekey: identifier of the node we wish to project onto the mesh
        #   plane_normal: the direction we wish to project the node onto the mesh
        #   cache: (optional) pass the cache to this function to enhance performance
      
        # Ouputs:
        #   new_coords: the new coordinates of the node now that it is on the mesh (or if we could not solve it will be the last iteration)
      
        # Adjust the solving parameters here:
        max_iterations = 20
        tol = .01
        plane_normal_offset = .3
        
       

        # Re-aquire openzinc variables
        fm = region.getFieldmodule()
        coordinates = fm.findFieldByName('coordinates')
        coordinates = coordinates.castFiniteElement()
        if cache == none:
          cache = fm.createFieldCache()

        # Create templates
        nodes = fm.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES)
        nodetemplate = nodes.createNodetemplate()
        nodetemplate.defineField(coordinates)
        nodetemplate.setValueNumberOfVersions(coordinates, -1, Node.VALUE_LABEL_VALUE, 1)

        # Create our first new node for the search
        old_node = nodes.findNodeByIdentifier(nodeKey)
        cache.setNode(old_node)
        [result, old_coords] = coordinates.evaluateReal(cache, 3)

        # Create our mesh search
        mesh = fm.findMeshByName('mesh3d')
        mesh_location = fm.createFieldStoredMeshLocation(mesh)
        found_mesh_location = fm.createFieldFindMeshLocation(coordinates, coordinates, mesh)
        found_mesh_location.setSearchMode(found_mesh_location.SEARCH_MODE_NEAREST)

        # shift the point in preparation for first iteration
        plane_norm = np.array(plane_normal)
        shifted_point = old_coords + plane_norm * plane_normal_offset
        coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, shifted_point.tolist())
        
        # initialise for our solver
        it = 1
        old_coords = shifted_point
        test_coords = [10, 10, 10]
        new_coords = shifted_point
        start3 = time.clock()
        
        
        while abs(np.linalg.norm(np.dot(test_coords, plane_norm) - np.dot(new_coords, plane_norm))) > tol:
            # ^^ test if x and y changes are within tolerence by testing the magnitude of the difference of the dot products
            #    in the direction of our normal

            # Find nearest mesh location using the evaluateMeshLocation function
            [el, coords] = found_mesh_location.evaluateMeshLocation(cache, 3)
            cache.setMeshLocation(el, coords)
            [result, mesh_coords] = coordinates.evaluateReal(cache, 3)

            # Update our search location by finding how far we have moved in the plane_normal directoin
            new_coords = old_coords + np.dot(mesh_coords - old_coords, plane_norm)*plane_norm
            cache.setNode(old_node)
            coordinates.setNodeParameters(cache, -1, Node.VALUE_LABEL_VALUE, 1, new_coords.tolist())

            # switch our test coordinates
            test_coords = old_coords
            old_coords = new_coords

            # Break in case we can not converge
            it += 1
            if it > max_iterations:
                print(f'Could not converge on node {nodeKey}')
                break
                
        print(f'Node {nodeKey} was solved in {it-1} iterations' )
        return(new_coords)