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'' (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.
hat the specific element needs to implement or specify are the following:
  * $\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 ''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>

Edit summary

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