This tutorial will describe the implementation of a standard plane stress triangle with linear approximation of the displacement field.

Each element, in the structural module (SM) needs to compute the following quantities:

Internal load vector $$\mathbf{f}_{\mathrm{int}} = \int_V \mathbf{B}^{\mathrm{T}} \mathbf{\sigma} \ \mathrm{d}V $$

External load vector $$\mathbf{f}_{\mathrm{ext}} = \int_V \mathbf{N}^{\mathrm{T}} \mathbf{b} \ \mathrm{d}V + \int_{\Gamma} \mathbf{N}^{\mathrm{T}} \mathbf{t} \ \mathrm{d}\Gamma $$

Tangent stiffness matrix $$\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.

This class implements the integration over the volume and the boundary. What the element needs to implement or specify is:

• \(\mathbf{N}\) and \(\mathbf{B}\) matrices
• An integration rule (number of integration points etc.)
• 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

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
}

Differential boundary element - computeEdgeVolumeAround

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();
}

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

this->interp.evaldNdx( dNdx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

$$\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 this matrix the \(\mathbf{B}\)-matrix can be constructed $$\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} $$

    // 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);

Auxilary metods that needs to be overloaded:

    /// 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 ); }