How to work with nodes and elements in python

When you work with python, in many applications you may need to access the nodes and elements data or you may want to apply some boundary conditions directly from python. In these cases you need to use some simple commands in your run file of python which we present them as below:

Contents 1 Nodes 1.1 Applying boundary conditions 2 Elements 2.1 Getting a value from integration points 2.2 Getting a value from integration points in C++

Nodes

Applying boundary conditions

To apply a boundary condition on some nodes, you should first detect your desirable nodes in the model and then save them in an array. We name this array as Dirirchlet_nodes and we first search all the nodes to detect some of them as our desirable nodes to apply a Dirichlet boundary condition. Dirichlet_nodes array will help us to access our desirable nodes faster, without searching all the nodes again. For example, we can add the following lines in our python run file, which we name it as Run.py:

Run.py:

Dirichlet_nodes = []  # Defining an array to store the nodes

for node in model_part.Nodes:

    if(node.Y > 100):
       node.Fix(DISPLACEMENT_X)
       node.Fix(DISPLACEMENT_Y)
       Dirichlet_nodes.append(node)

In this example we select the nodes whose Y coordinates are more than 100. To access the information of the nodes` coordinate, you can easily use the optional iterator node and then use the components of .X,.Y or .Z to access the X,Y and Z coordinates of each node respectively. model_part.Nodes is the interface to access the nodes information in python. For the selected nodes in this example, we want to fix the value of DISPLACEMENT_X and DISPLACEMENT_X with command .Fix and then we store each node in Dirichlet_nodes array with command .append. After fixing the nodes, we can apply any values for DISPLACEMENT_X and DISPLACEMENT_Y with the following commands:

for i in range(0,len(Dirichlet_nodes)):
    Dirichlet_nodes[i].SetSolutionStepValue(DISPLACEMENT_X,0,1) 
    Dirichlet_nodes[i].SetSolutionStepValue(DISPLACEMENT_Y,0,-1)

where we apply the values of 1 and -1 for DISPLACEMENT_X and DISPLACEMENT_Y respectively. In this way we just need to make some iterations through the loop with the length of len(Dirichlet_nodes) instead of the number of total nodes and it makes the implementation of the boundary condition faster specially when we have to apply this condition several times in our application, i.e. transient problems. To make the nodes free from previous boundary conditions we can use command .Free as

for i in range(0,len(Dirichlet_nodes)):
    Dirichlet_nodes[i].Free(DISPLACEMENT_X) 
    Dirichlet_nodes[i].Free(DISPLACEMENT_Y)

Elements

Getting a value from integration points

The typical example of working with integration points is to get stress or strain values from these points for an element. For the first step, you should get the current process information and put it in a variable such as proc_info:

proc_info = model_part.ProcessInfo 

then you can use the following loop to get a value such as Stress from integration point of each element in the model:

StressVector = [] # Defining an array to store stress components  

for element in model_part.Elements:
    StressVector = element.GetValuesOnIntegrationPoints(STRESS, proc_info)

where StressVector is the array to store stress components. The size of this vector depends on the output vector that you define in your element. With a 2D triangular linear element with one Gauss point, StressVector includes three components: two axial and one shear stresses. The important point is that when you work with an element with more than one Gauss point, the StressVector is not anymore a vector and it converts to a matrix with the number of rows equal to the number of Gauss points. The dimension of this matrix depends on your definition in the element and it can be changed. In order to clarify this definition, we mention one example where we define the value of stress in a triangular linear element with one Gauss point. The name of our element is Displacement2D and we present the part of C++ code inside this element with the filename of Displacement2D.cpp in custom_elements folder.

Getting a value from integration points in C++

For the first step, you need to add the definition name of STRESS in the function of GetValueOnIntegrationPoints() in your element:

Displacement2D.cpp:

void Displacement2D::GetValueOnIntegrationPoints(const Variable<Vector>& rVariable, std::vector<Vector>& rValues, const ProcessInfo& rCurrentProcessInfo)
{
 if(rVariable==STRESS)
    {
     CalculateOnIntegrationPoints( rVariable, rValues, rCurrentProcessInfo );
    }
}

The important output of this function is rValues which is defined as a Vector of vectors or a matrix. This output finally is transfered to python. For the next step, you should move on to the function of CalculateOnIntegrationPoints() and define a procedure to calculate stress components:

Displacement2D.cpp:

void Displacement2D::CalculateOnIntegrationPoints(const Variable<Vector>& rVariable, std::vector<Vector>& Output, const ProcessInfo& rCurrentProcessInfo)
{
  KRATOS_TRY;

  if(rVariable==STRESS)
  {

         array_1d<double,3> msStressVector;  // Defining a Stress Vector

         ## The code to calculate the msStressVector should be added here according to your application
         ##
         ##

         Temp.resize(msStressVector.size());
 
         for(unsigned int ii = 0; ii<msStressVector.size(); ii++)
            Temp[ii] = msStressVector[ii];

  Output[0] = Temp;
  }

  KRATOS_CATCH("")
}

To be continued ............