Tutorial:Creating the Conditions

From KratosWiki
Jump to: navigation, search

What you need to underestand of these lines

When the buider and solver is called in your strategy, it will loop through all the nodes and condition asking for their contribution to the system. This means that we need the following information:

- First of all, we need to know the contribution to the global system. In our case the 1x1 LHS matrix and the 1x1 RHS vector. Given the strategy we chose, to do we'll program this contribution in the CalculateLocalSystem. In this file we will perform all the necesary computations. We will return the rLeftHandSideMatrix (our LHS matrix) and the rRightHandSideVector (our RHS vector) . The rCurrentProcessInfo is an input used to gather general information, for example the delta_time in transcient problems. In this particular problem we are not using it due to its simplicity. SInce it is only a load, we will not modify the LHS matrix but only the RHS vector.

- Besides returning the 1x1 matrix and the RHS vector, we must tell the builder to which degrees of freedom each of the lines of this matrix must add its contribution. In our case it is pretty simple, the first line corresponds to the TEMPERATURE DoF of the first node, the seccond node to the TEMPERATURE dof of the second node and so on. This information is given to the builder and solver using both EquationIdVector and GetDofList.


In this particular case, since our condition is only linked to one Degree of freedom, it has a 1x1 matrix and a 1x1 RHS vector. If we were using a line condition (2 nodes), our matrix woud be 2x2 and the RHS 2x1.

In our case we will not use the CalculateRightHandSide, but the idea is basically the same as the CalculateLocalSystem


Creating a Custom Condition to incorporate the thermal loads: pointsource.h ; pointsource.cpp

The process of creating a custom condition is similar to creating an element, the easiest way is to copy the code from another condition and change it to suit our needs. Due to simplicity and similarities, the easiest files to start with are PointForce2D from the structural application. The difference between our recently created element and the condition we are about to create is that in this case we'll be loading data nodewisely.

POINT_HEAT_SOURCE is our variable containing the external thermal loads to our problem. To include them we'll create a new condition called pointsource . KRATOS will call the condition and assign it node by node to the system of equations. Since this is standard load, we only have to add the contribution to the right hand side (RHS) without any calculation. In pseudo code, our task is simply:

RHS (node) + = POINT_HEAT_SOURCE (node)

Still, note that the structre of a condition, as well as the one in elements, includes both the RHS and the left hand side matrix (LHS). This gives us flexibility and thus we can incorporate operations that also affect the matrix. We will create 2 files (pointsource.h and pointsource.cpp ), preferably inside the custom_conditions folder . Bear in mind that conditions are a class called by the solver , so you have to respect the structure and subroutines required by it. Again, the same rules that apply to elements. Again, in pseudocode:

class pointsource
           void CalculateLocalSystem
           void CalculateRightHandSide
           void EquationIdVector
           void GetDofList


Just as in the elements, to incorporate this condition we have to create it and register it in Kratos. Once registered, the solver will be the one calling it and ensambling the LHS and RHS, so there's no need to call it in the problemtype. As writen in the CMake file, the files below will have to be in the /custom_conditions folder. If it does not exist, create it and inside create 2 empty files. Below you can find both the .h and .cpp files:

pointsource.h

#if !defined(KRATOS_PointSource_CONDITION_H_INCLUDED )
#define  KRATOS_PointSource_CONDITION_H_INCLUDED

// External includes 
#include "boost/smart_ptr.hpp"

// Project includes
#include "includes/define.h"
#include "includes/serializer.h"
#include "includes/condition.h"
#include "includes/ublas_interface.h"
#include "includes/variables.h"

namespace Kratos
{
 class PointSource
  : public Condition
    {
    public:
      ///@name Type Definitions
      ///@{
      
       /// Counted pointer of PointForce2D
       KRATOS_CLASS_POINTER_DEFINITION(PointSource);
      
      /// Default constructor. 
      PointSource(IndexType NewId, GeometryType::Pointer pGeometry);

      PointSource(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties);

      /// Destructor.
      virtual ~PointSource();
      

      Condition::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes,  PropertiesType::Pointer pProperties) const;

      void CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo);
      
      void CalculateRightHandSide(VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo);

      //virtual void CalculateLeftHandSide(MatrixType& rLeftHandSideMatrix, ProcessInfo& rCurrentProcessInfo);
      
      void EquationIdVector(EquationIdVectorType& rResult, ProcessInfo& rCurrentProcessInfo);

      void GetDofList(DofsVectorType& ConditionalDofList,ProcessInfo& CurrentProcessInfo);

   protected:
 
 
   private:

	friend class Serializer;

	// A private default constructor necessary for serialization  
	PointSource() : Condition()
       {
       }
        

  }; // Class PointForce2D 

} //namespace kratos 
#endif


pointsource.cpp

// Project includes 
#include "includes/define.h"
#include "custom_conditions/pointsource.h"
#include "purediffusion_application.h"
#include "utilities/math_utils.h"

namespace Kratos
{
	//************************************************************************************
	//************************************************************************************
	PointSource::PointSource(IndexType NewId, GeometryType::Pointer pGeometry)
		: Condition(NewId, pGeometry)
	{		
		//DO NOT ADD DOFS HERE!!!
	}

	//************************************************************************************
	//************************************************************************************
	PointSource::PointSource(IndexType NewId, GeometryType::Pointer pGeometry,  PropertiesType::Pointer pProperties)
		: Condition(NewId, pGeometry, pProperties)
	{
	}
	Condition::Pointer PointSource::Create(IndexType NewId, NodesArrayType const& ThisNodes,  PropertiesType::Pointer pProperties) const
	{
		return Condition::Pointer(new PointSource(NewId, GetGeometry().Create(ThisNodes), pProperties));
	}

	PointSource::~PointSource()
	{
	}

	//************************************************************************************
	//************************************************************************************
	void PointSource::CalculateRightHandSide(VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
	{
		KRATOS_TRY
		if(rRightHandSideVector.size() != 1)
			rRightHandSideVector.resize(1,false);

		double load = GetGeometry()[0].GetSolutionStepValue(POINT_HEAT_SOURCE);
		rRightHandSideVector[0] = load;
		KRATOS_CATCH("")
	}

	//************************************************************************************
	//************************************************************************************
	void PointSource::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
	{
		KRATOS_TRY

		if(rLeftHandSideMatrix.size1() != 1)
			rLeftHandSideMatrix.resize(1,1,false);
		noalias(rLeftHandSideMatrix) = ZeroMatrix(1,1);
		if(rRightHandSideVector.size() != 1)
			rRightHandSideVector.resize(1,false);
		double load = GetGeometry()[0].GetSolutionStepValue(POINT_HEAT_SOURCE);
		rRightHandSideVector[0] = load;
		KRATOS_CATCH("")
	}


	//************************************************************************************
	//************************************************************************************
	void PointSource::EquationIdVector(EquationIdVectorType& rResult, ProcessInfo& CurrentProcessInfo)
	{
		int number_of_nodes = GetGeometry().PointsNumber();
		unsigned int index;
		unsigned int dim = 1;
		rResult.resize(number_of_nodes*dim);
		for (int i=0;i<number_of_nodes;i++)
		{
			index = i*dim;
			rResult[index] = (GetGeometry()[i].GetDof(TEMPERATURE).EquationId());			
		}
	}

	//************************************************************************************
	//************************************************************************************
	  void PointSource::GetDofList(DofsVectorType& ConditionalDofList,ProcessInfo& CurrentProcessInfo)
	{
		unsigned int dim = 1;
		ConditionalDofList.resize(GetGeometry().size()*dim);
		unsigned int index;
		for (unsigned int i=0;i<GetGeometry().size();i++)
		{
			
			index = i*dim;
			ConditionalDofList[index] = (GetGeometry()[i].pGetDof(TEMPERATURE));
		}
	}
} // Namespace Kratos




back to Kratos For Dummies

Personal tools
Categories