Skip to content

2D Triangle Element

Bases: Element

An up-to-3rd-order 2d triangular element with 3, 6 or 10 nodes respectively.

The node numbering follows the meshio conventions shown below:

1st order element::

2
|\
| \
|  \
|   \
|    \
0-----1

2nd order element::

2
|\
| \
5  4
|   \
|    \
0--3--1

3rd Order::

2
|\
| \
7  6
|   \
|    \
8  9  5
|      \
|       \
0--3--4--1
Source code in FEMpy/Elements/TriElement2D.py 
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
class TriElement2D(Element):
    """
    An up-to-3rd-order 2d triangular element with 3, 6 or 10 nodes respectively.

    The node numbering follows the meshio conventions shown below:

    1st order element::

        2
        |\\
        | \\
        |  \\
        |   \\
        |    \\
        0-----1

    2nd order element::

        2
        |\\
        | \\
        5  4
        |   \\
        |    \\
        0--3--1

    3rd Order::

        2
        |\\
        | \\
        7  6
        |   \\
        |    \\
        8  9  5
        |      \\
        |       \\
        0--3--4--1

    """

    def __init__(self, order=1, numStates=None, quadratureOrder=None):
        """Create a new 2d triangular finite element object

        Parameters
        ----------
        order : int, optional
            Element order, a first order quad has 4 nodes, 2nd order 9, 3rd order 16 etc, currently only orders 1-3 are
            supported, by default 1
        numStates : int, optional
            Number of states in the underlying PDE, by default 2
        quadratureOrder : int, optional
            Quadrature order to use for numerical integration, by default None, in which case a valid order for the
            chosen element order is used

        Raises
        ------
        ValueError
            Raises error if order is not 1, 2 or 3
        """
        if order not in [1, 2, 3]:
            raise ValueError("Triangular elements only support orders 1, 2 and 3")
        self.order = order
        numNodes = (order + 1) * (order + 2) // 2
        if quadratureOrder is None:
            shapeFuncOrder = order
            quadratureOrder = int(np.ceil((shapeFuncOrder + 1) / 2))

        super().__init__(numNodes, numDim=2, quadratureOrder=quadratureOrder, numStates=numStates)

        self.name = f"Order{self.order}-LagrangeTri"

        # --- Define parametric coordinate bounds ---
        self.paramCoordLowerBounds = -np.zeros(self.numDim)
        self.paramCoordLinearConstaintMat = np.array([1.0, 1.0])
        self.paramCoordLinearConstaintUpperBounds = 1.0
        self.paramCoordLinearConstaintLowerBounds = -np.inf

    def computeShapeFunctions(self, paramCoords):
        """Compute the shape function values at a given set of parametric coordinates

        Parameters
        ----------
        paramCoords : numPoint x numDim array
            Array of parametric point coordinates to evaluate shape functions at

        Returns
        -------
        N: numPoint x numNodes array
            Array of shape function values at the given parametric coordinates, N[i][j] is the value of the jth shape function at the ith parametric point
        """
        return LP.LagrangePolyTri(paramCoords[:, 0], paramCoords[:, 1], self.order)

    def computeShapeFunctionGradients(self, paramCoords):
        """Compute the derivatives of the shape functions with respect to the parametric coordinates at a given set of parametric coordinates

        Parameters
        ----------
        paramCoords : numPoint x numDim array
            Array of parametric point coordinates to evaluate shape function gradients at

        Returns
        -------
        NGrad: numPoint x numDim x numNodes array
            Shape function gradient values, NGrad[i][j][k] is the value of the kth shape function at the ith point w.r.t the kth
            parametric coordinate
        """
        return LP.LagrangePolyTriDeriv(paramCoords[:, 0], paramCoords[:, 1], self.order)

    def getIntegrationPointWeights(self, order=None):
        """Compute the integration point weights for a given quadrature order on this element

        Parameters
        ----------
        order : int
            Integration order

        Returns
        -------
        array of length numIntpoint
            Integration point weights
        """
        if order is None:
            order = self.quadratureOrder
        return getTriQuadWeights(order)

    def getIntegrationPointCoords(self, order=None):
        """Compute the integration point parameteric coordinates for a given quadrature order on this element

        Parameters
        ----------
        order : int
            Integration order

        Returns
        -------
        numIntpoint x numDim array
            Integration point coordinates
        """
        if order is None:
            order = self.quadratureOrder
        return getTriQuadPoints(order)

    def getReferenceElementCoordinates(self):
        """Get the node coordinates for the reference element, a.k.a the element on which the shape functions are defined

        Returns
        -------
        numNodes x numDim array
            Element node coordinates
        """
        if self.order == 1:
            return np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
        elif self.order == 2:
            return np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.5], [0.0, 0.5]])
        elif self.order == 3:
            return np.array(
                [
                    [0.0, 0.0],
                    [1.0, 0.0],
                    [0.0, 1.0],
                    [1 / 3, 0.0],
                    [2 / 3, 0.0],
                    [2 / 3, 1 / 3],
                    [1 / 3, 2 / 3],
                    [0.0, 2 / 3],
                    [0.0, 1 / 3],
                    [1 / 3, 1 / 3],
                ]
            )

    def getRandParamCoord(self, n, rng=None):
        """Generate a set of random parametric coordinates
        For a tri element we need u and v in range [0,1] and u + v <= 1, we can generate these points by generating
        random points in a square on the domain [0,1] and then reflecting any points outside the triangle to the inside.
        Parameters
        ----------
        n : int, optional
            number of points to generate, by default 1
        rng : numpy random Generator, optional
            Random number generator to use, useful for creating consistent test behaviour, by default None, in which
            case a new one is created for this call

        Returns
        -------
        paramCoords : n x numDim array
            isoparametric coordinates, one row for each point
        """
        if rng is None:
            rng = np.random.default_rng()
        coords = np.atleast_2d(rng.random((n, 2)))
        for i in range(n):
            coordSum = coords[i, 0] + coords[i, 1]
            if coordSum > 1:
                coords[i] = 1 - coords[i]
        return coords

__init__(order=1, numStates=None, quadratureOrder=None)

Create a new 2d triangular finite element object

Parameters

order : int, optional Element order, a first order quad has 4 nodes, 2nd order 9, 3rd order 16 etc, currently only orders 1-3 are supported, by default 1 numStates : int, optional Number of states in the underlying PDE, by default 2 quadratureOrder : int, optional Quadrature order to use for numerical integration, by default None, in which case a valid order for the chosen element order is used

Raises

ValueError Raises error if order is not 1, 2 or 3

Source code in FEMpy/Elements/TriElement2D.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
def __init__(self, order=1, numStates=None, quadratureOrder=None):
    """Create a new 2d triangular finite element object

    Parameters
    ----------
    order : int, optional
        Element order, a first order quad has 4 nodes, 2nd order 9, 3rd order 16 etc, currently only orders 1-3 are
        supported, by default 1
    numStates : int, optional
        Number of states in the underlying PDE, by default 2
    quadratureOrder : int, optional
        Quadrature order to use for numerical integration, by default None, in which case a valid order for the
        chosen element order is used

    Raises
    ------
    ValueError
        Raises error if order is not 1, 2 or 3
    """
    if order not in [1, 2, 3]:
        raise ValueError("Triangular elements only support orders 1, 2 and 3")
    self.order = order
    numNodes = (order + 1) * (order + 2) // 2
    if quadratureOrder is None:
        shapeFuncOrder = order
        quadratureOrder = int(np.ceil((shapeFuncOrder + 1) / 2))

    super().__init__(numNodes, numDim=2, quadratureOrder=quadratureOrder, numStates=numStates)

    self.name = f"Order{self.order}-LagrangeTri"

    # --- Define parametric coordinate bounds ---
    self.paramCoordLowerBounds = -np.zeros(self.numDim)
    self.paramCoordLinearConstaintMat = np.array([1.0, 1.0])
    self.paramCoordLinearConstaintUpperBounds = 1.0
    self.paramCoordLinearConstaintLowerBounds = -np.inf

computeShapeFunctionGradients(paramCoords)

Compute the derivatives of the shape functions with respect to the parametric coordinates at a given set of parametric coordinates

Parameters

paramCoords : numPoint x numDim array Array of parametric point coordinates to evaluate shape function gradients at

Returns

NGrad: numPoint x numDim x numNodes array Shape function gradient values, NGrad[i][j][k] is the value of the kth shape function at the ith point w.r.t the kth parametric coordinate

Source code in FEMpy/Elements/TriElement2D.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def computeShapeFunctionGradients(self, paramCoords):
    """Compute the derivatives of the shape functions with respect to the parametric coordinates at a given set of parametric coordinates

    Parameters
    ----------
    paramCoords : numPoint x numDim array
        Array of parametric point coordinates to evaluate shape function gradients at

    Returns
    -------
    NGrad: numPoint x numDim x numNodes array
        Shape function gradient values, NGrad[i][j][k] is the value of the kth shape function at the ith point w.r.t the kth
        parametric coordinate
    """
    return LP.LagrangePolyTriDeriv(paramCoords[:, 0], paramCoords[:, 1], self.order)

computeShapeFunctions(paramCoords)

Compute the shape function values at a given set of parametric coordinates

Parameters

paramCoords : numPoint x numDim array Array of parametric point coordinates to evaluate shape functions at

Returns

N: numPoint x numNodes array Array of shape function values at the given parametric coordinates, N[i][j] is the value of the jth shape function at the ith parametric point

Source code in FEMpy/Elements/TriElement2D.py 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def computeShapeFunctions(self, paramCoords):
    """Compute the shape function values at a given set of parametric coordinates

    Parameters
    ----------
    paramCoords : numPoint x numDim array
        Array of parametric point coordinates to evaluate shape functions at

    Returns
    -------
    N: numPoint x numNodes array
        Array of shape function values at the given parametric coordinates, N[i][j] is the value of the jth shape function at the ith parametric point
    """
    return LP.LagrangePolyTri(paramCoords[:, 0], paramCoords[:, 1], self.order)

getIntegrationPointCoords(order=None)

Compute the integration point parameteric coordinates for a given quadrature order on this element

Parameters

order : int Integration order

Returns

numIntpoint x numDim array Integration point coordinates

Source code in FEMpy/Elements/TriElement2D.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
def getIntegrationPointCoords(self, order=None):
    """Compute the integration point parameteric coordinates for a given quadrature order on this element

    Parameters
    ----------
    order : int
        Integration order

    Returns
    -------
    numIntpoint x numDim array
        Integration point coordinates
    """
    if order is None:
        order = self.quadratureOrder
    return getTriQuadPoints(order)

getIntegrationPointWeights(order=None)

Compute the integration point weights for a given quadrature order on this element

Parameters

order : int Integration order

Returns

array of length numIntpoint Integration point weights

Source code in FEMpy/Elements/TriElement2D.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def getIntegrationPointWeights(self, order=None):
    """Compute the integration point weights for a given quadrature order on this element

    Parameters
    ----------
    order : int
        Integration order

    Returns
    -------
    array of length numIntpoint
        Integration point weights
    """
    if order is None:
        order = self.quadratureOrder
    return getTriQuadWeights(order)

getRandParamCoord(n, rng=None)

Generate a set of random parametric coordinates For a tri element we need u and v in range [0,1] and u + v <= 1, we can generate these points by generating random points in a square on the domain [0,1] and then reflecting any points outside the triangle to the inside. Parameters

n : int, optional number of points to generate, by default 1 rng : numpy random Generator, optional Random number generator to use, useful for creating consistent test behaviour, by default None, in which case a new one is created for this call

Returns

paramCoords : n x numDim array isoparametric coordinates, one row for each point

Source code in FEMpy/Elements/TriElement2D.py 
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
def getRandParamCoord(self, n, rng=None):
    """Generate a set of random parametric coordinates
    For a tri element we need u and v in range [0,1] and u + v <= 1, we can generate these points by generating
    random points in a square on the domain [0,1] and then reflecting any points outside the triangle to the inside.
    Parameters
    ----------
    n : int, optional
        number of points to generate, by default 1
    rng : numpy random Generator, optional
        Random number generator to use, useful for creating consistent test behaviour, by default None, in which
        case a new one is created for this call

    Returns
    -------
    paramCoords : n x numDim array
        isoparametric coordinates, one row for each point
    """
    if rng is None:
        rng = np.random.default_rng()
    coords = np.atleast_2d(rng.random((n, 2)))
    for i in range(n):
        coordSum = coords[i, 0] + coords[i, 1]
        if coordSum > 1:
            coords[i] = 1 - coords[i]
    return coords

getReferenceElementCoordinates()

Get the node coordinates for the reference element, a.k.a the element on which the shape functions are defined

Returns

numNodes x numDim array Element node coordinates

Source code in FEMpy/Elements/TriElement2D.py 
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
def getReferenceElementCoordinates(self):
    """Get the node coordinates for the reference element, a.k.a the element on which the shape functions are defined

    Returns
    -------
    numNodes x numDim array
        Element node coordinates
    """
    if self.order == 1:
        return np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
    elif self.order == 2:
        return np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [0.5, 0.0], [0.5, 0.5], [0.0, 0.5]])
    elif self.order == 3:
        return np.array(
            [
                [0.0, 0.0],
                [1.0, 0.0],
                [0.0, 1.0],
                [1 / 3, 0.0],
                [2 / 3, 0.0],
                [2 / 3, 1 / 3],
                [1 / 3, 2 / 3],
                [0.0, 2 / 3],
                [0.0, 1 / 3],
                [1 / 3, 1 / 3],
            ]
        )