oofem-wiki

You've loaded an old revision of the document! If you save it, you will create a new version with this data.

====== 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.

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

<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 is evaluated using the elements' interpolator
<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 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 \\
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} \\
\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>

<code c++>
void
BasicElement :: computeBmatrixAt(GaussPoint *gp, FloatMatrix &answer, int li, int ui)
{
    /* Compute the [3x6] strain-displacement matrix {B} for the element,
     * evaluated at the given gp.
     * {B}*{a} should provide the strains {eps} in Voigt form 
     * {eps} = {eps_xx, eps_yy, gam_xy}^T with {a} being the 
     * solution vector of the element.
     */
    
    /* Evaluate the derivatives of the shape functions at the position of the gp.     
     * dNdx = [dN1/dx1 dN1/dx2
     *         dN2/dx1 dN2/dx2
     *         dN3/dx1 dN3/dx2]
    */  
    FloatMatrix dNdx; 
    this->interp.evaldNdx( dNdx, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );

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

Part of header file...
<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>

Edit summary

Note: By editing this page you agree to license your content under the following license: CC Attribution-Share Alike 4.0 International