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

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 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} $$

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

BasicElement.C

BasicElement.h

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.