\(
   \renewcommand{\b}[1]{\mathbf{#1}}
\)
$$\b{a}_{\rm{i}} \rm{b} b$$

====== Implementation of a linear plane stress triangle ======
This tutorial will describe the implementation of a standard plane stress triangle with linear approximation of the displacement field.
For quasi-static problems the following set of equations needs to be solved
$$\mathbf{f}_{\mathrm{int}} = \mathbf{f}_{\mathrm{ext}}$$
with the internal and external force vectors respectively

$$\b{f}_{\rm{int}} = \int_V \b{B}^{\rm{T}} \b{\sigma} \ \rm{d}V $$


$$\b{f}_{\rm{ext}} = \int_V \b{N}^{\rm{T}} \b{b} \ \rm{d}V 
+ \int_{\Gamma} \b{N}^{\rm{T}} \b{t} \ \rm{d}\Gamma $$

The default solution procedure for solving the equations are a Newton-Rapshon scheme and for this the tangent stiffness matrix is needed,
computed as 
$$\mathbf{K} = \int_V \mathbf{B}^{\mathrm{T}} \mathbf{D} \mathbf{B} \ \mathrm{d}V $$

Since the FE equations for all standard continuum elements (2D plane stress/strain, 3D, etc.) are equal, much 
of the element implementations are placed in the base class ''NonLinStructuralElement'' (or ''StructuralElement'' for purely linear problems).
For example, this base class implements methods to compute the internal and external load vector together with the tangent stiffnes matrix (all as defined above).
The main question that arises is then: what must a new element implement? In short it is everything that is element specific, such as
the following:
  * \(\mathbf{N}\) and \(\mathbf{B}\) matrices
    * These depends on the number of nodes and the number o dofs stored in each node.
  * An integration rule
    * This defines what type of integration alogorithm (Gauss, Newton-Cotes, Lobatto), the number of integration points, 
    * if there should be several algorithms to support reduced integration for example.
  * Compute the differential volume element \(\Delta V \) and \(\Delta \Gamma \)

In addition to the above several auxilary methods needs to be implemented such as:
  * Number of nodes and dofs
  * What type of dofs that exists in the nodes
  * Read specific element parameters from the input file (not standard ones...)

All these topics will be considered in this tutorial.


Differential volume element - ''computeVolumeAround''
Differential boundary element

<code c++>
double 
BasicElement :: computeVolumeAround(GaussPoint *gp)
{
    // Computes the volume element dV associated with the given gp.

    double weight = gp->giveWeight();
    FloatArray *lCoords = gp->giveNaturalCoordinates(); // local/natural coords of the gp (parent domain)
    double detJ = fabs( this->interp.giveTransformationJacobian( *lCoords, FEIElementGeometryWrapper(this) ) );
    double thickness = this->giveCrossSection()->give(CS_Thickness, gp); // the cross section keeps track of the thickness
    
    return detJ * thickness * weight; // dV
}
</code>

Differential boundary element - ''computeEdgeVolumeAround''
<code c++>
double
BasicElement :: computeEdgeVolumeAround(GaussPoint *gp, int iEdge)
{
    /* Returns the line element ds associated with the given gp on the specific edge.
     * Note: The name is misleading since there is no volume to speak of in this case. 
     * The returned value is used for integration of a line integral (external forces).
     */
    double detJ = this->interp.edgeGiveTransformationJacobian( iEdge, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
    return detJ * gp->giveWeight();
}
</code>


Computation of the \(\mathbf{B}\)-matrix:
This method is called ''computeBmatrixAt''

The derivatives of the shape functions ''dNdx'' are evaluated using the elements' interpolator class
<code c++>
this->interp.evaldNdx( dNdx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
</code>

$$\frac{d \mathbf{N}}{d \mathbf{x}}
= \begin{pmatrix} 
\frac{\mathrm{d}N_1}{\mathrm{d}x_1} & \frac{\mathrm{d}N_1}{\mathrm{d}x_2} \\
\frac{\mathrm{d}N_2}{\mathrm{d}x_1} & \frac{\mathrm{d}N_2}{\mathrm{d}x_2} \\
\frac{\mathrm{d}N_3}{\mathrm{d}x_1} & \frac{\mathrm{d}N_3}{\mathrm{d}x_2} 
\end{pmatrix}
$$ 

From these derivatives the \(\mathbf{B}\)-matrix can be constructed which for a three noded triangle with two dofs per node looks like
$$\mathbf{B}
= \begin{pmatrix} 
\frac{\mathrm{d}N_1}{\mathrm{d}x_1} & 0 & \frac{\mathrm{d}N_2}{\mathrm{d}x_1} & 0 & \frac{\mathrm{d}N_3}{\mathrm{d}x_1} & 0 \\
\frac{\mathrm{d}N_1}{\mathrm{d}x_1} & 0 & \frac{\mathrm{d}N_2}{\mathrm{d}x_1} & 0 & \frac{\mathrm{d}N_3}{\mathrm{d}x_1} & 0\\
\frac{\mathrm{d}N_1}{\mathrm{d}x_2} & \frac{\mathrm{d}N_1}{\mathrm{d}x_1} & \frac{\mathrm{d}N_2}{\mathrm{d}x_2} 
& \frac{\mathrm{d}N_2}{\mathrm{d}x_1} & \frac{\mathrm{d}N_3}{\mathrm{d}x_2} & \frac{\mathrm{d}N_3}{\mathrm{d}x_1} \\
\end{pmatrix}
$$ 
<code c++>
    // Construct the B-matrix
    answer.resize(3, 6);

    answer.at(1, 1) = dNdx.at(1, 1);
    answer.at(1, 3) = dNdx.at(2, 1);
    answer.at(1, 5) = dNdx.at(3, 1);

    answer.at(2, 2) = dNdx.at(1, 2);
    answer.at(2, 4) = dNdx.at(2, 2);
    answer.at(2, 6) = dNdx.at(3, 2);

    answer.at(3, 1) = dNdx.at(1, 2);
    answer.at(3, 2) = dNdx.at(1, 1);
    answer.at(3, 3) = dNdx.at(2, 2);
    answer.at(3, 4) = dNdx.at(2, 1);
    answer.at(3, 5) = dNdx.at(3, 2);
    answer.at(3, 6) = dNdx.at(3, 1);
</code>





Auxilary metods that needs to be overloaded:
<code c++>
    
    /// Constructor
    BasicElement(int n, Domain * d);    
    /// Destructor.
    virtual ~BasicElement() { }         

    virtual FEInterpolation *giveInterpolation() const;
    virtual int computeNumberOfDofs() { return 6; }
    virtual void giveDofManDofIDMask(int inode, IntArray &) const;
    virtual double computeVolumeAround(GaussPoint *gp);

    
    virtual const char *giveInputRecordName() const { return _IFT_BasicElement_Name; }
    virtual const char *giveClassName() const { return "BasicElement"; }
    virtual IRResultType initializeFrom(InputRecord *ir);
    virtual MaterialMode giveMaterialMode() { return _PlaneStress; }
    virtual int testElementExtension(ElementExtension ext) { return ( ( ext == Element_EdgeLoadSupport ) ? 1 : 0 ); }
</code>

<file - BasicElement.C>
</file>

<file - BasicElement.h>
</file>


===== Problem representation - Engineering model =====

The concept "engineering model" is an abstraction for the problem under
consideration. It represents the type of analysis to be performed (e.g. static structural, transient heat flow, etc.).
The base class ''EngngModel'' declares and implements the basic services for assembling
characteristic components and services for starting the solution step and
its termination. Derived classes ``know'' the form of governing
equation and the physical meaning of  particular components. 
They are responsible for forming the governing equation for each solution
step,  usually by summing contributions from particular elements and
nodes.ecific load type dependent
services and implement all necessary services.