# Magnetostatic Application Elements

## Introduction

The magnetostatic Kratos element is the implementation code of the Elemental Finite Element Formulation for the Magnetostatic PDE. This element can be used for any already existing application in Kratos, by including it in their corresponding files.

In this point we will need just to know the name of the application for their corresponding files ("kMagnetostatic"), and the name of the application methods "MagnetostaticApplication".

Our first element in Kratos will be the implementation of the Poisson's equation in Magnetostatic, so we will call it "magnetostatic_2d". The class name for the element will be called "Magnetostatic2D".

The most basic magnetostatic element implemented in Kratos only needs two files:

• magnetostatic_2d.h
• magnetostatic_2d.cpp

A quick way to create an element is to select one of the existing ones close to our formulation and modifying some specific parts. Magnetostatic Kratos element is almost the same to the electrostatic one, so we strongly recommend you to use this as a reference template. We are going to show the code of these two files with some comments in those lines you probably would want to customise.

The following code includes the most usual sections needed to operate with a Kratos element. Most of the code is common in all the elements and it doesn't need any change. Specific customisation is remarked using bold characters, therefore in general terms you should pay attention just to these parts and keep the rest unchanged.

First lines: legal issues

The first lines are just the typical legal issues related to Kratos as a framework to create simulation programs. You should keep them as they are. The only sentence to customise is the name of the application (in bold characters) in the third line, KratosR1MagnetostaticApplication:

/*
==============================================================================
KratosR1MagnetostaticApplication
A library based on:
Kratos
A General Purpose Software for Multi-Physics Finite Element Analysis
Version 1.0 (Released on march 05, 2007).

Pooyan Dadvand, Riccardo Rossi, Javier Mora
pooyan@cimne.upc.edu
rrossi@cimne.upc.edu
mora@cimne.upc.edu
- CIMNE (International Center for Numerical Methods in Engineering),
Gran Capita' s/n, 08034 Barcelona, Spain

Permission is hereby granted, free  of charge, to any person obtaining
a  copy  of this  software  and  associated  documentation files  (the
"Software"), to  deal in  the Software without  restriction, including
without limitation  the rights to  use, copy, modify,  merge, publish,
distribute,  sublicense and/or  sell copies  of the  Software,  and to
permit persons to whom the Software  is furnished to do so, subject to
the following condition:

Distribution of this code for  any  commercial purpose  is permissible
ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNERS.

The  above  copyright  notice  and  this permission  notice  shall  be
included in all copies or substantial portions of the Software.

THE  SOFTWARE IS  PROVIDED  "AS  IS", WITHOUT  WARRANTY  OF ANY  KIND,
EXPRESS OR  IMPLIED, INCLUDING  BUT NOT LIMITED  TO THE  WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT  SHALL THE AUTHORS OR COPYRIGHT HOLDERS  BE LIABLE FOR ANY
CLAIM, DAMAGES OR  OTHER LIABILITY, WHETHER IN AN  ACTION OF CONTRACT,
TORT  OR OTHERWISE, ARISING  FROM, OUT  OF OR  IN CONNECTION  WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

==============================================================================
*/


The following lines (comments again) indicate the author's name (your name), date and hour of the last revision, and revision number, as usual in other codes.

//
//   Project Name:        Kratos
//   Last modified by:    $Author: jmora$
//   Date:                $Date: 2010-02-02 11:58:38$
//   Revision:            $Revision: 1.0$
//
//


In order to avoid class redefinition, you have to define a label for your element definition, as follows:

#if !defined(KRATOS_MAGNETOSTATIC_2D_ELEM_H_INCLUDED )
#define  KRATOS_MAGNETOSTATIC_2D_ELEM_H_INCLUDED


the #endif is at the end of the code (the rest of the lines must be included between these ones):

#endif // KRATOS_MAGNETOSTATIC_2D_ELEM_H_INCLUDED  defined


Now you have to add the "include" files (copy as it is):

// System includes

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

// Project includes
#include "includes/define.h"
#include "includes/element.h"
#include "includes/ublas_interface.h"
#include "includes/variables.h"


Declare the namespace Kratos (copy at it is):

namespace Kratos
{

///@name Kratos Globals
///@{

///@}
///@name Type Definitions
///@{

///@}
///@name  Enum's
///@{

///@}
///@name  Functions
///@{

///@}
///@name Kratos Classes
///@{

/// Short class definition.
/** Detail class definition.
*/


The namespace is closed at the end of the file (just before the #endif):

}  // namespace Kratos.


Our class definition starts now. We will call it "Magnetostatic2D":

 class Magnetostatic2D
: public Element
{
public:
///@name Type Definitions
///@{

/// Counted pointer of Magnetostatic2D
KRATOS_CLASS_POINTER_DEFINITION(Magnetostatic2D);

///@}
///@name Life Cycle
///@{


Constructor and destructor declaration (you should just change the class name):

    /// Default constructor.
Magnetostatic2D(IndexType NewId, GeometryType::Pointer pGeometry);
Magnetostatic2D(IndexType NewId, GeometryType::Pointer pGeometry,  PropertiesType::Pointer pProperties);

/// Destructor.
virtual ~ Magnetostatic2D();

///@}
///@name Operators
///@{

///@}
///@name Operations
///@{


And we add the class methods we need to perform the implementation of the element. By now we will just declare the methods, leaving a more detailed explanation in the implementation file (.cpp). Note that for this simple element you don't need change anything:

      Element::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& ElementalDofList,ProcessInfo& CurrentProcessInfo);

void InitializeSolutionStep(ProcessInfo& CurrentProcessInfo);

void CalculateOnIntegrationPoints( const Variable<array_1d<double,3> >& rVariable,
std::vector<array_1d<double,3> >& rValues,
const ProcessInfo& rCurrentProcessInfo);

void GetValueOnIntegrationPoints(const Variable<array_1d<double,3>>& rVariable, std::vector<array_1d<double,3> >& rValues, const ProcessInfo& rCurrentProcessInfo);



Additional lines not important for this element (you don't need to modify anything):

    ///@}
///@name Access
///@{

///@}
///@name Inquiry
///@{

///@}
///@name Input and output
///@{

///@}
///@name Friends
///@{

///@}

protected:
///@name Protected static Member Variables
///@{

///@}
///@name Protected member Variables
///@{

///@}
///@name Protected Operators
///@{

///@}
///@name Protected Operations
///@{

///@}
///@name Protected  Access
///@{

///@}
///@name Protected Inquiry
///@{

///@}
///@name Protected LifeCycle
///@{

///@}


The declaration of the class's attributes we will know to perform the element computation (note again the use of the class name):

   private:
///@name Static Member Variables
///@{
static boost::numeric::ublas::bounded_matrix<double,3,2> msDN_DX;
static boost::numeric::ublas::bounded_matrix<double,2,2> msD;
static array_1d<double,3> msN; //dimension = number of nodes
static array_1d<double,3> ms_temp; //dimension = number of nodes
static array_1d<double,3> point_sources; //dimension = number of nodes


and to close the class definition:

     ///@}
///@name Member Variables
///@{

///@}
///@name Private Operators
///@{

///@}
///@name Private Operations
///@{

///@}
///@name Private  Access
///@{

///@}
///@name Private Inquiry
///@{

///@}
///@name Un accessible methods
///@{

/// Assignment operator.
// Magnetostatic2D& operator=(const Poisson2D& rOther);

/// Copy constructor.
// Magnetostatic2D(const Poisson2D& rOther);

///@}

}; // Class Magnetostatic2D


To finalise the header code, we should add (no changes are needed):

 ///@}

///@name Type Definitions
///@{

///@}
///@name Input and output
///@{

/// input stream function

/// output stream function

///@}


## Element implementation methods file: magnetostatic_2d.cpp

First lines: legal issues

The first lines are just again the typical legal issues related to Kratos as a framework to create simulation programs. You should keep them as they are. The only sentence to customise is the name of the application (in bold characters) in the third line, KratosR1MagnetostaticApplication:

/*
==============================================================================
KratosR1MagnetostaticApplication
A library based on:
Kratos
A General Purpose Software for Multi-Physics Finite Element Analysis
Version 1.0 (Released on march 05, 2007).

Pooyan Dadvand, Riccardo Rossi, Javier Mora
pooyan@cimne.upc.edu
rrossi@cimne.upc.edu
mora@cimne.upc.edu
- CIMNE (International Center for Numerical Methods in Engineering),
Gran Capita' s/n, 08034 Barcelona, Spain

Permission is hereby granted, free  of charge, to any person obtaining
a  copy  of this  software  and  associated  documentation files  (the
"Software"), to  deal in  the Software without  restriction, including
without limitation  the rights to  use, copy, modify,  merge, publish,
distribute,  sublicense and/or  sell copies  of the  Software,  and to
permit persons to whom the Software  is furnished to do so, subject to
the following condition:

Distribution of this code for  any  commercial purpose  is permissible
ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNERS.

The  above  copyright  notice  and  this permission  notice  shall  be
included in all copies or substantial portions of the Software.

THE  SOFTWARE IS  PROVIDED  "AS  IS", WITHOUT  WARRANTY  OF ANY  KIND,
EXPRESS OR  IMPLIED, INCLUDING  BUT NOT LIMITED  TO THE  WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT  SHALL THE AUTHORS OR COPYRIGHT HOLDERS  BE LIABLE FOR ANY
CLAIM, DAMAGES OR  OTHER LIABILITY, WHETHER IN AN  ACTION OF CONTRACT,
TORT  OR OTHERWISE, ARISING  FROM, OUT  OF OR  IN CONNECTION  WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

==============================================================================
*/


The following lines (comments again) indicate the author's name (your name), date and hour of the last revision, and revision number, as usual in other codes.

//
//   Project Name:        Kratos
//   Last modified by:    $Author: jmora$
//   Date:                $Date: 2010-02-02 11:58:38$
//   Revision:            $Revision: 1.0$
//
//


Now you have to add the "include" files. Note that you have to change the element header file name (magnetostatic_2d.h) and the application header file name where you want to include your element (in this case, kMagnetostatic.h, not yet created). The related subdirectories will be explained in the application section:

// System includes

// External includes

// Project includes

#include "includes/define.h"
#include "custom_elements/magnetostatic_2d.h"
#include "kMagnetostatic.h"
#include "utilities/math_utils.h"
#include "utilities/geometry_utilities.h"


Declare the namespace Kratos (copy at it is). List the static variables, and create the main class methods as follows. Just note that you have to be sure that you are using your element class name (Magnetostatic2D in this example, in bold characters):

namespace Kratos
{
namespace Magnetostatic2Dauxiliaries
{
//static variables
boost::numeric::ublas::bounded_matrix<double,3,2> msDN_DX = ZeroMatrix(3,2);
boost::numeric::ublas::bounded_matrix<double,2,3> msB = ZeroMatrix(2,3);
boost::numeric::ublas::bounded_matrix<double,2,2> msD = ZeroMatrix(2,2);

array_1d<double,3> msN = ZeroVector(3); //dimension = number of nodes
array_1d<double,3> ms_temp = ZeroVector(3); //dimension = number of nodes
array_1d<double,3> point_sources = ZeroVector(3); //dimension = number of nodes

array_1d<double,3> surface_sources = ZeroVector(3); //dimension = number of nodes
}
using  namespace Magnetostatic2Dauxiliaries;

//************************************************************************************
//************************************************************************************
Magnetostatic2D::Magnetostatic2D(IndexType NewId, GeometryType::Pointer pGeometry)
: Element(NewId, pGeometry)
{
}

//************************************************************************************
//************************************************************************************
Magnetostatic2D::Magnetostatic2D(IndexType NewId, GeometryType::Pointer pGeometry,  PropertiesType::Pointer pProperties)
: Element(NewId, pGeometry, pProperties)
{

}

Element::Pointer Magnetostatic2D::Create(IndexType NewId, NodesArrayType const& ThisNodes,  PropertiesType::Pointer pProperties) const
{
return Element::Pointer(new Magnetostatic2D(NewId, GetGeometry().Create(ThisNodes), pProperties));
}

Magnetostatic2D::~Magnetostatic2D()
{
}


Remember to close the namespace at the end of the file.

}  // namespace Kratos.


Now, it's the critical moment to transfer your Magnetostatic Finite Element Formulation to the Kratos code. Some preliminary comments are:

• Kratos allows different ways to implement the formulation, therefore don't worry if you find different Kratos elements with different implementation approaches, even if you leave some methods unused, or parts without contents. Just check the coherence of your formulation.
• Kratos has been designed to allow you to focus your attention to the element formulation. That is, again, don't worry about the other programming aspects of your simulation program, such as the matrix assembly, solvers, etc... Of course, you are able to customise them, but for basic applications is not necessary at all.

Copy the following code as it is, with two exceptions (in bold characters):

• use your element class name (Magnetostatic2D in this example);
• your variable names (we will use MAGNETOSTATIC_POTENTIAL for the unknown, MAGNETIC_PERMEABILITY for the material property, MAGNETIC_FIELD_INTENSITY and MAGNETIC_FLUX_DENSITY for the gradient of the unknown, INFINIT_COEFFICIENT for the infinit boundary condition, MAGNETOSTATIC_SURFACE_CURRENT and MAGNETOSTATIC_POINT_CURRENT for the source terms);
	//************************************************************************************
//************************************************************************************
void Magnetostatic2D::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
KRATOS_TRY

const unsigned int number_of_points = GetGeometry().size();

if(rLeftHandSideMatrix.size1() != number_of_points)
rLeftHandSideMatrix.resize(number_of_points,number_of_points,false);

if(rRightHandSideVector.size() != number_of_points)
rRightHandSideVector.resize(number_of_points,false);

//getting data for the given geometry
double Area;
GeometryUtils::CalculateGeometryData(GetGeometry(), msDN_DX, msN, Area);

array_1d<double,3> permeability = GetProperties()[MAGNETIC_PERMEABILITY];
msD(0,0)=1/permeability[1];
msD(1,1)=1/permeability[0];
msD(1,0)=0.0;
msD(0,1)=0.0;

//double surface_sources = (this)->GetValue(MAGNETOSTATIC_SURFACE_CURRENT);
surface_sources[0] = (this)->GetValue(MAGNETOSTATIC_SURFACE_CURRENT);
surface_sources[1] = (this)->GetValue(MAGNETOSTATIC_SURFACE_CURRENT);
surface_sources[2] = (this)->GetValue(MAGNETOSTATIC_SURFACE_CURRENT);

// main loop
const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints();

noalias(rLeftHandSideMatrix) = prod(msDN_DX,Matrix(prod(msD,trans(msDN_DX))));


Note that this line is related to the following term of the stiffness matrix:

$\mathbf{K}^{(e)}= \int_{\Omega^{(e)}} \mathbf{B^T} \mathbf{\varepsilon} \mathbf{B} \partial \Omega^{(e)}$

		rLeftHandSideMatrix *= Area;

noalias(rRightHandSideVector) = surface_sources*Area/3.0;


and here we add the surface term to the right side of the expression:

$\mathbf{f}^{(e)}= \int_{\Omega^{(e)}} \mathbf{N^T} \rho_v \partial \Omega^{(e)}$
		//subtracting the dirichlet term
// RHS -= LHS*MAGNETOSTATIC_POTENTIALs
for(unsigned int iii = 0; iii<number_of_points; iii++)
ms_temp[iii] = GetGeometry()[iii].FastGetSolutionStepValue(MAGNETOSTATIC_POTENTIAL);
noalias(rRightHandSideVector) -= prod(rLeftHandSideMatrix,ms_temp);

KRATOS_CATCH("");
}



The code is self-explanatory, but just a few comments to link it with our basic finite element formulation:

• CalculateGeometryData obtains msDN_DX, msN and the elemental Area, that is, the B and N matrix;
• msD is the material properties matrix (the D matrix);
• With prod(msDN_DX,Matrix(prod(msD,trans(msDN_DX)))) we obtain the stiffness matrix (the K matrix);
• With noalias(rRightHandSideVector) = surface_sources we add the source_source contribution to the source vector f;
• Finally, with prod(rLeftHandSideMatrix,ms_temp) we compute the Dirichlet contribution to the source vector f;

The following code is not used for this specific element (keep at it is, just changing the class name):

	//************************************************************************************
//************************************************************************************
void Magnetostatic2D::CalculateRightHandSide(VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
KRATOS_ERROR(std::logic_error,  "method not implemented" , "");
}

//************************************************************************************
//************************************************************************************
// this subroutine calculates the nodal contributions for the explicit steps of the
// fractional step procedure
void Magnetostatic2D::InitializeSolutionStep(ProcessInfo& CurrentProcessInfo)
{
KRATOS_TRY

KRATOS_CATCH("");
}


The next lines give us the unknows (you have just change the class and variables' names):

	//************************************************************************************
//************************************************************************************
void Magnetostatic2D::EquationIdVector(EquationIdVectorType& rResult, ProcessInfo& CurrentProcessInfo)
{
unsigned int number_of_nodes = GetGeometry().PointsNumber();
if(rResult.size() != number_of_nodes)
rResult.resize(number_of_nodes,false);

for (unsigned int i=0;i<number_of_nodes;i++)
rResult[i] = GetGeometry()[i].GetDof(MAGNETOSTATIC_POTENTIAL).EquationId();
}

//************************************************************************************
//************************************************************************************
void Magnetostatic2D::GetDofList(DofsVectorType& ElementalDofList,ProcessInfo& CurrentProcessInfo)
{
unsigned int number_of_nodes = GetGeometry().PointsNumber();
if(ElementalDofList.size() != number_of_nodes)
ElementalDofList.resize(number_of_nodes);

for (unsigned int i=0;i<number_of_nodes;i++)
ElementalDofList[i] = GetGeometry()[i].pGetDof(MAGNETOSTATIC_POTENTIAL);

}


And, finally, the following methods compute the solution for the gradient of the unknowns:

	//************************************************************************************
//************************************************************************************
void Magnetostatic2D::CalculateOnIntegrationPoints(const Variable<array_1d<double,3>>& rVariable, std::vector<array_1d<double,3> >& Output, const ProcessInfo& rCurrentProcessInfo)
{
IntegrationMethod mThisIntegrationMethod;

mThisIntegrationMethod= GetGeometry().GetDefaultIntegrationMethod();
const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints(mThisIntegrationMethod);
Vector vect_tmp(3);
vect_tmp[0] = GetGeometry()[0].FastGetSolutionStepValue(MAGNETOSTATIC_POTENTIAL);
vect_tmp[1] = GetGeometry()[1].FastGetSolutionStepValue(MAGNETOSTATIC_POTENTIAL);
vect_tmp[2] = GetGeometry()[2].FastGetSolutionStepValue(MAGNETOSTATIC_POTENTIAL);

array_1d<double,3> permeability = GetProperties()[MAGNETIC_PERMEABILITY];
msD(0,0)=1/permeability[1];
msD(1,1)=1/permeability[0];
msD(1,0)=0.0;
msD(0,1)=0.0;

double Area;
GeometryUtils::CalculateGeometryData(GetGeometry(), msDN_DX, msN, Area);

if(Output.size() != GetGeometry().IntegrationPoints(mThisIntegrationMethod).size())
Output.resize(GetGeometry().IntegrationPoints(mThisIntegrationMethod).size());

for(unsigned int PointNumber = 0; PointNumber<integration_points.size(); PointNumber++)
{
if(rVariable==MAGNETIC_FLUX_DENSITY)
{
//KRATOS_WATCH("GiD Post Magnetostatic - calc - magn - field");
noalias(Output[PointNumber])=-prod(trans(msDN_DX),vect_tmp);
Output[PointNumber][2]=0.0;
//				KRATOS_WATCH(Output[PointNumber]);
}


CHECK THIS!!!! It's wrong! That is:

$B_x = -\frac{\partial A_z(x,y)}{\partial y} \qquad B_y = -\frac{\partial A_z(x,y)}{\partial x}$

			else if(rVariable==MAGNETIC_FIELD_INTENSITY)
{
//KRATOS_WATCH(rValues[PointNumber]);
noalias(Output[PointNumber])=-prod(prod(msD,trans(msDN_DX)),vect_tmp);
//				KRATOS_WATCH(Output[PointNumber]);
}


equivalent to:

$H_x = -\frac{1}{\mu_x} \frac{\partial A_z(x,y)}{\partial y} \quad H_y = -\frac{1}{\mu_y} \frac{\partial A_z(x,y)}{\partial x}$
			else
{
//				KRATOS_WATCH("GiD Post Magnetostatic - calc - else");
Output[PointNumber][0]=msD(0,0);
Output[PointNumber][1]=msD(1,1);
Output[PointNumber][2]=0.0;
//				KRATOS_WATCH(Output[PointNumber]);
}
}

}

//************************************************************************************
//************************************************************************************

void Magnetostatic2D::GetValueOnIntegrationPoints(const Variable<array_1d<double,3> >& rVariable,
std::vector<array_1d<double,3> >& rValues, const ProcessInfo& rCurrentProcessInfo)
{
//KRATOS_WATCH("GiD Post Magnetostatic - GetValueOnIntegrationPoints");

if(rVariable==MAGNETIC_FLUX_DENSITY)
{
//			KRATOS_WATCH("GiD Post Magnetostatic - get - magn-flux");
CalculateOnIntegrationPoints( rVariable, rValues, rCurrentProcessInfo );
}
else if(rVariable==MAGNETIC_FIELD_INTENSITY)
{
//			KRATOS_WATCH("GiD Post Magnetostatic - get - magn-inten");
CalculateOnIntegrationPoints( rVariable, rValues, rCurrentProcessInfo );
}
else
{
//			KRATOS_WATCH("GiD Post Magnetostatic - get - else");
CalculateOnIntegrationPoints( rVariable, rValues, rCurrentProcessInfo );
}
}