User:Eml5526.s11.team7/HW6

Problem Description
Similar to HW5.{1,3,7}, solve G1DM1.0/D1b using Quadratic Lagrangian Element Basis Function (QLEBF) with uniform discretization (equidistant element nodes) until convergence of $${{u}^{h}}(0.5)$$ to $$O({{10}^{-6}})$$. ( $$nel=2,4,6,8,...$$).


 * 1. For $$nel=2$$, compute $$\mathbf{\tilde{K}}=\sum\limits_{e=1}^{2}$$, with $${{\mathbf{\tilde{K}}}^{e}}$$ by (1)p.30-5, display $${{\mathbf{\tilde{K}}}^{e}},e=1,2$$.


 * 2. Compute $${{\mathbf{k}}^{e}}$$, $${{\mathbf{L}}^{e}}$$, for $$e=1,2$$.


 * 3. Compute $${{\tilde{\Kappa }}^{e}}={{\mathbf{L}}^{e}}^{T}{{\mathbf{k}}^{e}}{{\mathbf{L}}^{e}}$$, for $$e=1,2$$, compare to the result by (1).


 * 4. Plot all QLEBF for $$nel=3$$.


 * 5. Plot $$u_^{h}$$ vs. $$u$$, $$\left[ u_^{h}(0.5)-u(0.5) \right]$$ vs. $$\tilde{n}$$

Part 1: Compute K tilde Matrix: $$\displaystyle \underline{\tilde K}=\sum_{e=1}^2\underline{\tilde K^e}$$
For element 1, the 3 basis functions are:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_1^1 = \frac \cdot \frac = 8{x^2} - 6x + 1

$$     (6.1.1)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_2^1 = \frac \cdot \frac = - 16{x^2} + 8x

$$     (6.1.2)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_3^1 = \frac \cdot \frac = 8{x^2} - 2x

$$     (6.1.3) By (1)p.30-5,
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{ij}^e = \int\limits_ {b{{_i^e}^\prime }(x){a_2}(x)b{{_j^e}^\prime }(x)dx}

$$     (6.1.4) Thus, we have
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{11}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( {16x - 6} \right)dx} = \frac{6}

$$     (6.1.5)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{12}^1 = K_{21}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( { - 32x + 8} \right)dx} =  - \frac{3}

$$     (6.1.6)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{13}^1 = K_{31}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} = \frac{6}

$$     (6.1.7)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{22}^1 = \int\limits_0^{0.5} {\left( { - 32x + 8} \right)\left( {2 + 3x} \right)\left( { - 32x + 8} \right)dx} = \frac{3}

$$     (6.1.8)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{23}^1 = K_{32}^1 = \int\limits_0^{0.5} {\left( { - 32x + 8} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} =  - \frac{3}

$$     (6.1.9)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{33}^1 = \int\limits_0^{0.5} {\left( {16x - 2} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} = \frac{6}

$$     (6.1.10) Then
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde K}}^1} = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{89/6}&0&0 \\   0&0&0&0&0 \\   0&0&0&0&0 \end{array}} \right]

$$     (6.1.11) Similar to element 1, when  e=2
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde K}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0 \\  0&0&0&0&0 \\   0&0&{107/6}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.1.12) Thus,
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{\mathbf{\tilde K}} = \sum\limits_{e = 1}^2  = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{98/3}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.1.13)
 * 
 * }

Part 2: Compute element k and L matrices: $$\displaystyle \underline{k}^e,\underline{L}^e$$, for e=1,2
By (1)p.31-3,
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

k_{ij}^e = \int\limits_ {b{{_i^e}^\prime }(x){a_2}(x)b{{_j^e}^\prime }(x)dx}

$$     (6.1.14) For element 1
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{k}}^1} = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6} \\  { - 38/3}&{88/3}&{ - 50/3} \\   {11/6}&{ - 50/3}&{89/6} \end{array}} \right]

$$     (6.1.15)
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 1&0&0&0&0 \\  0&1&0&0&0 \\   0&0&1&0&0 \end{array}} \right]

$$     (6.1.16) For element 2
 * 
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{k}}^2} = \left[ {\begin{array}{ccccccccccccccc} {107/6}&{ - 62/3}&{17/6} \\  { - 62/3}&{136/3}&{ - 74/3} \\   {17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.1.17)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{L}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&1&0&0 \\  0&0&0&1&0 \\   0&0&0&0&1 \end{array}} \right]

$$     (6.1.18)
 * <p style="text-align:right">
 * }

Part 3: Compute element K tilde matrix from k tilde and scatter matrices: $$\displaystyle \underline{\tilde K}^e=\underline{L}^{e^T}\underline{k}^e\underline{L}^e$$ for e=1,2
For element 1
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde {\rm K}}}^1} = {{\mathbf{L}}^1}^T{{\mathbf{k}}^1}{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 1&0&0 \\  0&1&0 \\   0&0&1 \\   0&0&0 \\   0&0&0 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6} \\  { - 38/3}&{88/3}&{ - 50/3} \\   {11/6}&{ - 50/3}&{89/6} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 1&0&0&0&0 \\  0&1&0&0&0 \\   0&0&1&0&0 \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{89/6}&0&0 \\   0&0&0&0&0 \\   0&0&0&0&0 \end{array}} \right]

$$     (6.1.19) For element 2
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde {\rm K}}}^2} = {{\mathbf{L}}^2}^T{{\mathbf{k}}^2}{{\mathbf{L}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0 \\  0&0&0 \\   1&0&0 \\   0&1&0 \\   0&0&1 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {107/6}&{ - 62/3}&{17/6} \\  { - 62/3}&{136/3}&{ - 74/3} \\   {17/6}&{ - 74/3}&{131/6} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0&0&1&0&0 \\  0&0&0&1&0 \\   0&0&0&0&1 \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0 \\  0&0&0&0&0 \\   0&0&{107/6}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.1.20) The result is exactly same as achieved in part 1.
 * <p style="text-align:right">
 * }

Part 4: Plot all QLEBF for nel=3
The plot of QLEBF for nel = 3:

Note: Parts 1-4 in 6.2 are exactly the same as 1-4 in 6.1. Jiang has done the work to solve 1-4 in part 6.2.

Part 5: Plot u approximate vs u exact, and the error in u approximate at .5 vs number of elements: $$\displaystyle u^hvsu$$, $$\displaystyle [u^h(.5)-u(.5)]vs\tilde n$$
The exact solution is:

$$u = - \frac{5}{x^2} + \frac{5}{9}x + \frac\log \left( {3x + 2} \right) + 4 - \frac\log \left( 2 \right)$$

The plots of $$u_^{h}$$ vs. $$u$$(nel=14)



The plots of $$\left[ u_^{h}(0.5)-u(0.5) \right]$$ vs. $$\tilde{n}$$

Hua

Problem Description
Similar to HW5.{2,4,8}, solve G1DM1.0/D1b using Quadratic Lagrangian Element Basis Function (QLEBF) with uniform discretization (equidistant element nodes) until convergence of $${{u}^{h}}(0.5)$$ to $$O({{10}^{-6}})$$. ( $$nel=2,4,6,8,...$$).


 * 1. For $$nel=2$$, compute $$\mathbf{\tilde{K}}=\sum\limits_{e=1}^{2}$$, with $${{\mathbf{\tilde{K}}}^{e}}$$ by (1)p.30-5, display $${{\mathbf{\tilde{K}}}^{e}},e=1,2$$.


 * 2. Compute $${{\mathbf{k}}^{e}}$$, $${{\mathbf{L}}^{e}}$$, for $$e=1,2$$.


 * 3. Compute $${{\tilde{\Kappa }}^{e}}={{\mathbf{L}}^{e}}^{T}{{\mathbf{k}}^{e}}{{\mathbf{L}}^{e}}$$, for $$e=1,2$$, compare to the result by (1).


 * 4. Plot all QLEBF for $$nel=3$$.


 * 5. Plot $$u_^{h}$$ vs. $$u$$, $$\left[ u_^{h}(0.5)-u(0.5) \right]$$ vs. $$\tilde{n}$$

Part 1: Compute K tilde Matrix: $$\displaystyle \underline{\tilde K}=\sum_{e=1}^2\underline{\tilde K^e}$$
For element 1, the 3 basis functions are:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_1^1 = \frac \cdot \frac = 8{x^2} - 6x + 1

$$     (6.2.1)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_2^1 = \frac \cdot \frac = - 16{x^2} + 8x

$$     (6.2.2)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

b_3^1 = \frac \cdot \frac = 8{x^2} - 2x

$$     (6.2.3) By (1)p.30-5,
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{ij}^e = \int\limits_ {b{{_i^e}^\prime }(x){a_2}(x)b{{_j^e}^\prime }(x)dx}

$$     (6.2.4) Thus, we have
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{11}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( {16x - 6} \right)dx} = \frac{6}

$$     (6.2.5)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{12}^1 = K_{21}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( { - 32x + 8} \right)dx} =  - \frac{3}

$$     (6.2.6)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{13}^1 = K_{31}^1 = \int\limits_0^{0.5} {\left( {16x - 6} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} = \frac{6}

$$     (6.2.7)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{22}^1 = \int\limits_0^{0.5} {\left( { - 32x + 8} \right)\left( {2 + 3x} \right)\left( { - 32x + 8} \right)dx} = \frac{3}

$$     (6.2.8)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{23}^1 = K_{32}^1 = \int\limits_0^{0.5} {\left( { - 32x + 8} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} =  - \frac{3}

$$     (6.2.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

K_{33}^1 = \int\limits_0^{0.5} {\left( {16x - 2} \right)\left( {2 + 3x} \right)\left( {16x - 2} \right)dx} = \frac{6}

$$     (6.2.10) Then
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde K}}^1} = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{89/6}&0&0 \\   0&0&0&0&0 \\   0&0&0&0&0 \end{array}} \right]

$$     (6.2.11) Similar to element 1, when  e=2
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde K}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0 \\  0&0&0&0&0 \\   0&0&{107/6}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.12) Thus,
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{\mathbf{\tilde K}} = \sum\limits_{e = 1}^2  = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{98/3}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.13)
 * <p style="text-align:right">
 * }

Part 2: Compute element k and L matrices: $$\displaystyle \underline{k}^e,\underline{L}^e$$, for e=1,2
By (1)p.31-3,
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

k_{ij}^e = \int\limits_ {b{{_i^e}^\prime }(x){a_2}(x)b{{_j^e}^\prime }(x)dx}

$$     (6.2.14) For element 1
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{k}}^1} = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6} \\  { - 38/3}&{88/3}&{ - 50/3} \\   {11/6}&{ - 50/3}&{89/6} \end{array}} \right]

$$     (6.2.15)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 1&0&0&0&0 \\  0&1&0&0&0 \\   0&0&1&0&0 \end{array}} \right]

$$     (6.2.16) For element 2
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{k}}^2} = \left[ {\begin{array}{ccccccccccccccc} {107/6}&{ - 62/3}&{17/6} \\  { - 62/3}&{136/3}&{ - 74/3} \\   {17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.17)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{L}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&1&0&0 \\  0&0&0&1&0 \\   0&0&0&0&1 \end{array}} \right]

$$     (6.2.18)
 * <p style="text-align:right">
 * }

Part 3: Compute element K tilde matrix from k tilde and scatter matrices: $$\displaystyle \underline{\tilde K}^e=\underline{L}^{e^T}\underline{k}^e\underline{L}^e$$ for e=1,2
For element 1
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde {\rm K}}}^1} = {{\mathbf{L}}^1}^T{{\mathbf{k}}^1}{{\mathbf{L}}^1} = \left[ {\begin{array}{ccccccccccccccc} 1&0&0 \\  0&1&0 \\   0&0&1 \\   0&0&0 \\   0&0&0 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6} \\  { - 38/3}&{88/3}&{ - 50/3} \\   {11/6}&{ - 50/3}&{89/6} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 1&0&0&0&0 \\  0&1&0&0&0 \\   0&0&1&0&0 \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} {65/6}&{ - 38/3}&{11/6}&0&0 \\  { - 38/3}&{88/3}&{ - 50/3}&0&0 \\   {11/6}&{ - 50/3}&{89/6}&0&0 \\   0&0&0&0&0 \\   0&0&0&0&0 \end{array}} \right]

$$     (6.2.19) For element 2
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{{\mathbf{\tilde {\rm K}}}^2} = {{\mathbf{L}}^2}^T{{\mathbf{k}}^2}{{\mathbf{L}}^2} = \left[ {\begin{array}{ccccccccccccccc} 0&0&0 \\  0&0&0 \\   1&0&0 \\   0&1&0 \\   0&0&1 \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} {107/6}&{ - 62/3}&{17/6} \\  { - 62/3}&{136/3}&{ - 74/3} \\   {17/6}&{ - 74/3}&{131/6} \end{array}} \right]\left[ {\begin{array}{ccccccccccccccc} 0&0&1&0&0 \\  0&0&0&1&0 \\   0&0&0&0&1 \end{array}} \right] = \left[ {\begin{array}{ccccccccccccccc} 0&0&0&0&0 \\  0&0&0&0&0 \\   0&0&{107/6}&{ - 62/3}&{17/6} \\   0&0&{ - 62/3}&{136/3}&{ - 74/3} \\   0&0&{17/6}&{ - 74/3}&{131/6} \end{array}} \right]

$$     (6.2.20) The result is exactly same as achieved in part 1.
 * <p style="text-align:right">
 * }

Part 4: Plot all QLEBF for nel=3
The plot of QLEBF for nel = 3:

Part 5: Plot u approximate vs u exact, and the error in u approximate at .5 vs number of elements: $$\displaystyle u^hvsu$$, $$\displaystyle [u^h(.5)-u(.5)]vs\tilde n$$
The exact solution is:

$$u = - \frac{5}{x^2} + \frac{5}{9}x - \frac\log \left( {3x + 2} \right) + \frac{5}{12} + \frac\log \left( 5 \right)+4-\frac{5}{9}$$

The plots of $$u_^{h}$$ vs. $$u$$,



The plots of $$\left[ u_^{h}(0.5)-u(0.5) \right]$$ vs. $$\tilde{n}$$

Jiang

Part1: Building Geometry and Meshing
After installing CGX and CCX launch BConverge and save a blank file with the name beam.fbd:



By saving the file with this extension the file will automatically format itself (similar to when you save a .m file or a .f file in an IDE). To get CGX running with this as the input go to Tools->Pre Process, giving the following screen:



Now enter in the following information into beam.fbd:



Geometry can be created through the command line in an interactive format, similar to running from the command line in Matlab, or geometry can be created from a batch file, similar to a Matlab script. The batch file is the easiest way to assign geometry where using the command line to adjust certain parameters thereafter is convenient. With this batch file we’ve created two points, p1 and p2, and one line, l1, which connects the two points. As before, run this batch file, Tools->Pre Process. You’ll notice that the graph has changed:



One can now play around with the various functions the mouse has inside the window where the line is plotted. Holding down the left mouse button and moving the mouse allows rotation of the object (in this case the line). Holding down the middle button (pressing in on the scroll wheel on most mouses) and moving the mouse allows for zooming in and out. Lastly, holding the right mouse button down and moving the mouse allows for translational movement of the object. Clicking outside of the black bordered in object window allows for adjustment to what is being viewed on screen. These different options (seen when one presses the left mouse button outside of the black boarder in the Calculix Graphix (CGX) window) can also be executed through the command line. In fact, there are even more adjustments available through the command line than the drop down menu in the window.

Now the batch file will be adjusted for a slightly less elementary plot:



Run the file (Tools->PreProcess, note may have to exit out of the plotted window if already open).



Now, work will be done through the command line while being in the CGX window. By clicking inside the black boarder of the CGX window, type “plot pa all” and “plus la all”, then return. The following will be seen:



“Plot pa all” plots all the points in the set named all (or all the points) with their names. The “plus la all” command plots with the points all the lines belonging to the set all (so all the lines) and their names. If the command typed would have just been “plot la all” after “plot pa all”, the points would have been removed and just the lines would have been added. The plus allows for both to be plotted at once.

“Plot pa all” plots all the points in the set named all (or all the points) with their names. The “plus la all” command plots with the points all the lines belonging to the set all (so all the lines) and their names. If the command typed would have just been “plot la all” after “plot pa all”, the points would have been removed and just the lines would have been added. The plus allows for both to be plotted at once.

To further develop our FEA abilities with Calculix, were going to create a surface. This is done using the sweep command. To do this alter the batch file to be the following:



Clicking inside the black boarder of the CGX window (so as to allow for commands to be entered), enter the following commands, “plot pa all”, “plus la all”, “plus sa all”, giving the following:



So the seta command has created a set of lines named lines to sweep. The swep command sweeps the lines through a translational movement in the y direction (0,10,0) with the last number creating 10 divisions and storing these new lines and points into the set labeled sweplines. This correspondingly has created surfaces and lines due to the sweeping. To see the lines with their divisions, type the command “plot ld all” (just as we have before, clicking inside the black boarder on the Calculix window). This will give you:



The 25 and 10 are the corresponding divisions, so we now have a mesh with a resolution of 100x10. To now use what we’ve learned to make a 3D body, sweep in the z direction. This is done with the following batch file:



Running Pre processing, and adjusting the background color through Viewing->Toggle Background Color on the drop down menu outside the black boarder of the CGX window gives:



A body has now been created and we can start meshing! This will be done by issuing the following commands, “elty all he8”, “mesh all”, “plot m all”. The following will be seen:



Note the viewing effects have been turned off in the same way they were turned on. What has just happened is the elty all he8 command told CGX it’ll be meshing the set all with he8 elements. There are a variety of elements available to mesh with on CGX (be2:2node 1d linear beam element to pe 15:15 node 3D quadratic penta element). More complicated meshing can be done through another program like gmsh or netgen and then import the mesh to add boundary conditions and loads.

Now we will distinguish between model edges and element edges. To remove model edges use the drop down menu on the CGX window Viewing->Toggle Model Edges. To see the element edges clearly, use the drop down menu Viewing->DOTS, then Viewing ->Toggle Element Edges. You’ll get:



Or if you zoom in using the middle mouse button (clicking and holding down the scroll wheel and moving the mouse) you’ll see:



This shows that we’ve successfully created the mesh!

Part 2: Exporting Mesh,Loads, and Boundary Conditions
Calculix FEA Beam Part 2: Exporting Mesh, Loads, and Boundary Conditions We’ve now created a model of a beam. This model is made up of elements defined as he8 in Calculix, and has a mesh volume of 100x10x1. This next part of the series will create input files for CCX to solve based off of the created beam and the loads and boundary conditions we will apply.

Reopen the beam.fbd file if it isn’t already. We will be creating input files in Abaqus format, which is how CCX reads the input. Note, CCX borrows Abaqus syntax, but is built from scratch. Alter the beam.fbd file according to the following picture:



What has happened to the file is we’ve used the command set open (seto) to open a set named beam. Then everything we have done to build the beam in the first part of the tutorial is included, and then the set is closed through the set close command (setc). Note, a new set could be opened and the same descriptions done to create different components of a model. The lines following setc beam are used to code the type of element to be used and consequently build the mesh. Then send beam abq exports the created mesh from the previous commands into the Abaqus format. This creates a file in the same directory that is being worked in with the name beam.msh. The last two lines of beam.fbd are used to change the display window (CGX window). The first rotates the view to the z direction and the frame command allows the entire object to be seen in the frame. Using Tools->Pre Process will create the beam.msh file and give the following:



If one was to open beam.msh he or she would see the corresponding values the ccx reads as the mesh created for the figure shown. One could enter these explicitly in a .inp file which CCX would read, orthis file can just be loaded into the .inp file using the *INCLUDE command. A note about the distinction between CCX and CGX commands. CCX commands start with an * and are followed with an all uppercase word. All lower case letters and no * is a CGX command.

Now we will create the loads and boundary conditions that will also be used by CCX to produce our results. This will be done through a combination of using the command line and the CGX window. Pulling up the CGX window and clicking inside the border in order to issue commands as done previously, enter the following commands: “plot n all”, “rot –z”. You’ll see the following:



To remove the black border around our model in the CGX window, use the same steps previously used (Viewing->Toggle Model Edges). To have the image be less crowded use Viewing->DOTS. What the previously entered plot n all command achieved was to plot the nodes (n) for the set all (so all nodes). Then rot –z command has the same result as before (rotating view to be at z axis). Toggle back on the model edges and move to the left hand portion of the beam (Viewing->Toggle Model Edges; click and hold right mouse button and move mouse to view left hand portion of beam; click and hold middle mouse button (scroll wheel) and move mouse to zoom in) giving the following:



We will now be specifying this edge’s nodes to be fixed. Previously seta was used as a command based way to add items to a set. Now, the command qadd will allow for a graphical way to select items to be added to a set. This allows one to not need to know the explicit values of an item in order to add it to a set. So now to add the required items, the command qadd fixed can be issued (after clicking inside the black border in the CGX window). This will allow addition of items until q is pressed on the keyboard. Any item inside the black border at the tip of the mouse’s arrow will be selected as seen below. This can be altered in order to create a larger selection window (add more items at one time). To do this, after entering “qadd fixed” for the command press “r” making the current mouse position the top left of the future selection rectangle; move the mouse pointer to what will be the bottom right hand corner and press “r” again. A new selection rectangle has just been defined. Since we are now going to add a bulk amount of nodes, we must be in mode:a. This is done by entering “a” on the keyboard while still in qadd mode (after entering qadd fixed command). Move the new rectangle over the nodes to be added, as seen here:



With the proper nodes selected, enter n. If one returns to the command window he/she can see that the nodes have been added. These nodes have been added to the set fixed. Now, return to the CGX window to enter “q” in order to leave qadd.

If someone would like to confirm that he/she added the correct nodes, this can be done through the following commands: “view edge off”; “plus l beam w”; “plus n fixed g”. These commands first turn off the model edges, then plus l beam adds the lines from the set beam in white (“w”) and plus n fixed g adds the nodes (“n”) from the set fixed in green (“g”). After entering these commands one will see (after toggling background color as before):



What we have now done is select a set of nodes to assign a boundary condition to. For this set we will be fixing the nodes. To do this enter the following command (after clicking inside the black border of CGX window) “send fixed abq spc 123”. What this does is take the set fixed and (“fixed”) and uses single point constraint (“spc”) on each degree of freedom (“123”). Nodes have only three degrees of translational freedom, but for elements (things that can rotate etc) have 6 degrees of freedom. Execution of this send command has created the file for the boundary conditions fixed_123.bou. Now a load will be added to the bar. To view element faces enter the command “plot f beam”. Toggling on the element edges will also allow an easier view of the element faces. Move the beam to see the top part of it (right and left mouse buttons accordingly):



Use qadd again but this time the command is “qadd load”, then select the face accordingly as seen below and press “f” for face to add to the load set:



Now press q to exit qadd. Similarly to creating the boundary condition file, we will now create a load file. This is done through the following command “send load abq pres 10” where send creates the load.dlo file, load specifies the set, the format is Abaqus (“abq”) the ouput is a pressure (“pres”) with magnitude 10. As a check that the the Abaqus files have bee n created check the folder of your path:



Note specifying the load on this face could have been done through the command line as well. In order to do that one would need to the element number, and the node numbers. The commands to give the same results but with output file loadnodes.frc are “seta loadnodes n 1676 1673”; “send loadnodes abq force 0 -5 0”.

We have now created the necessary files to run CCX!

Part 3: Writing an Input File for the CCX Solver
Calculix FEA Beam Part 3: Writing an Input File for the CCX Solver Using the information from Part 1 and Part 2 we will now build an input file for the CCX solver. Opening a new file titled beam.inp:



Enter in the following information:



So we’ve created an input file for our beam that is fixed and one end and loaded at the other (cantilever beam). Lines 1 through 9 build the attributes for the model we’ll be running calculations on; lines 10 through 18 are the loading step portions. It is possible to print out the results to a file through the loading step. Note, these two parts make up the components of the input file, the model and the step.

Additional physical information is provided in lines 6 through 9. First the material is assigned the name EL, and then the elastic properties (through the *ELASTIC command) are loaded. Last, the command *SOLID SECTION assigns the material (we gave the name EL to) to the element set Ebeam. Looking now at the step portion of the file, *STEP begins this portion and *STATIC tells the CCX that we are working with a quasi-static case. So, the *STEP portion allows for calculation of now steady state cases where one can “march” forward in time. Here though, the *STATIC command allows for ignoring of the mass inertia of the structure. The command *INCLUDE includes the load file previously created to include the applied force at the end of the beam. The *NODE FILE and *EL FILE create files for the output of the step. Now one just needs to process this file. This is done through Tools->Solve. Upon successful completion one will see:
 * HEADING allows for a line of comments underneath it. *INCLUDE command allows for that previously collected data of the model to be loaded into the input file (note could explicitly type the information here, or using the include command to copy and paste the data to this input file).  Note the beam.msh file which is loaded for information about the model is already in the correct format for CCX to read.  However, now we will load the boundary condition file (fixed_123.bou) which is not already in the required format for CCX to read.  This is taken into account through the command *BOUNDARY preceding the loading of the boundary condition file.



A file titled beam.frb has been created in the working directory. We can post process this file through use of cgx. Go to Tools->Post Process. The following will pop up:



Using the dropdown menu inside the CGX window go to Datasets->2 STRESS 1.000000. Then go to Datasets->Entity->Mises. To see the Von Mises stress distribution at step one:



Part 4: Post Processing Results in CGX
In addition to creating the model and mesh to be used by CCX, CGX can provide post calculation processing. Open beam.inp created previously and go to Tools->Post Processing in order to open CGX with the file beam.frd. The following window will pop up:



Left clicking outside the black border in the CGX window to pull up the dropdown menu, one will notice there is now a section entitled Datasets. This option is only available when working on a .frb file with simulation results. Our simulation has two possible fields for visualization, displacement and stress. To see these things go to the dropdown menu in the CGX window, Datasets->Entity->All :



Since the displacement in this simulation is very small, it won’t be visible on the figure unless we make it so. First, to add displacement go to the dropdown menu outside the black border box in the CGX window, Viewing->Toggle Add-Dsisplcement. This is done through scaling. To scale up our displacement, enter the command “scal d 10000” giving:



To return back to regular scaling the command is simply “scal d”. To return to a known state and have the default gray beam, go to the dropdown menu outside the black border box in the CGX window, then Viewing->Show All Elements with Light.

A movie of the beam can also be recorded. First, by playing with the various features under Animate in the CGX dropdown menu one can achieve the animation they wish to display. Then the command “movi delay 0.01” is used to set the time between images taken of the beam movement. Depending on how many steps one has chosen by going to the drop down menu in the CGX window, Animate->Steps per Period-># Steps has selected, issue the command “movi start” and once the output reaches the number of chosen steps, issue the command “movi stop”. Then, if one chose 12 steps, issue the command “movi make 1 12”. This will create a movie.gif file. The command “movi clean” will delete the additional frames taken from the hard drive. The command prompt results will show the following:



Next a cross section of the beam can be examined to see the effects of the applied conditions. This is done through the cut command on the drop down menu obtained by clicking outside the black boarder in the CGX window. First, one can return to the original scaling through the command “scal d”. Then by going to Cut->Node 1 the model will be split up into its elements. Focus then on the plane you with to form with cut running through the bar. Cut requires selection of three nodes which will determine the place “cutting” through the model. So, after zooming in (middle mouse button (scroll wheel) and dragging) on the desired three nodes, left click on the first one. Then go to Cut->Node 2 and select the second node, and lastly Cut->Node 3 to select the third node. After this you’ll see:



Further analysis can be done with the use of the graph command. So multiple graphs can be made we’ll make a set that will be stored as a sequence. The storing the set as a sequence preserves the order in which things were added. This is important for example when plotting displacement vs x as will be done. When adding things to a sequence through qadd using the bulk addition, the order is based off of the element number, which will then plot the values in this order. This would not create a plot in order on x then. So, to create our plot of displacement vs x, first ensure the proper data is with the figure, Dataset->Entity->All. Then to better see the nodes, Viewing->DOTS, Viewing->Toggle Element Edges, “rot –z” to view from z axis, and “rot u 45” to allow for the top line to be at x=0,y=10, and z=0. Next, issue the following command “qadd graphseq s”. The s at the end indicates that qadd is adding a sequence. This will allow the graphing afterwards to plot correctly. Resize the square if needed through the “r” command (“r” once to set upper left of selection square, move cursor, “r” again to set lower right of selection square) then add the desired nodes (at the top) one at a time by placing the selection box over the node and pressing “n”. Once the nodes on screen have been selected, press “q” to exit out of qadd, use the right hand mouse button to translate to the new nodes to be added, and repeat (“qadd graphseq s”, “r” resize “r”, “n”). Once all the data has been added to the grapseq, enter the following command “graph graphseq l” where l is for length which has graph calculate the distance between successive nodes for the plot. The following plot will be made:



Lastly, to see the Values in the CGX window, enter the command “plot nv all” to see:



Text book (Fish and Belytschko) Problem 6.4.1

 * {| style="width:100%" border="0"

$$ Given\; the\; vector\; field\; {\color{blue}{\mathbf{ q_x=-y^2, q_y=-2xy}}} \;on\; the\; domain\; shown\; figure.\; Verify\; the\; divergence\; theorem. $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$Divergence\; Theorem\; is\; :\;
 * style="width:10%; |
 * style="width:10%; |

\color{blue}{\mathbf{{\int_{\Omega}\overrightarrow{\nabla}\cdot \overrightarrow{q}d\Omega} = {\oint_{\Gamma} \overrightarrow{q} \cdot \overrightarrow{n} d\Gamma}}} $$     (6.4.1)
 * <p style="text-align:right">
 * }

Solution

 * {| style="width:100%" border="0"

$$ \overrightarrow{\nabla}\cdot\overrightarrow{q}=\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}=0+\left(-2x\right)=-2x $$     (6.4.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \int_{\Omega}\overrightarrow{\nabla}\cdot \overrightarrow{q}d\Omega=\int_{-1}^{1}{\left( \int_{-1}^{1}-2xy\cdot dy \right) \cdot dx} = \int_{-1}^{1}\left(-xy^{2}\right)_{-1}^{1}dx=\int_{-1}^{1} 0 \cdot dx = \color{blue}\mathbf{0} $$     (6.4.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \overrightarrow{q}=q_x\overrightarrow{i}+q_y\overrightarrow{j} $$     (6.4.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \overrightarrow{q}=-y^2 \overrightarrow{i}-2xy\overrightarrow{j} $$     (6.4.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \oint{\overrightarrow{q}\cdot \overrightarrow{n}d\Gamma}={\int_{AB}2xydx} + {\int_{BC}-y^{2}dy}+{\int_{CD}2xydx}+{\int_{DA}y^{2}dy} $$     (6.4.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \oint{\overrightarrow{q}\cdot \overrightarrow{n}d\Gamma}=\left( x^{2}y\right)_{-1}^{1}+\left(\frac{-y^{3}}{3} \right)_{-1}^{1}+\left( x^{2}y\right)_{1}^{-1}+\left( \frac{y^{3}}{3}\right)_{-1}^{1}=\color{blue}\mathbf{0} $$     (6.4.7)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$\color{blue}6.4.7==6.4.3 \;Thus\; divergence\; theorem\; is\; verified.$$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Text book (Fish and Belytschko) Problem 6.4.2

 * {| style="width:100%" border="0"

$$ Given\; the\; vector\; field\; {\color{blue}{\mathbf{ q_x=3x^{2}y+y^3, q_y=3x+y^{3}}}} \;on\; the\; domain\; shown\; in\; figure.\; Verify\; the\; divergence\; theorem. $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \;The\; curved\; boundary\; of\; the\; domain\; is\; a\; parabola. $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Solution

 * {| style="width:100%" border="0"

$$ \overrightarrow{\nabla}\cdot\overrightarrow{q}=\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}=6xy+3y^{2} $$     (6.4.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$Equation\; of\; the\; given\; parabola\; y=4-x^2$$
 * style="width:10%; |
 * style="width:10%; |

(6.4.9)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$\int_{\Omega}\overrightarrow{\nabla}\cdot \overrightarrow{q}d\Omega=$$ $$\int_{-2}^{2}{\left( \int_{0}^{4-x^{2}}\left( 6xy+3y^{2}\right)\cdot dy \right) \cdot dx}=$$ $$\int_{-2}^{2}3x\left( 4-x^{2}\right)^{2}+$$$$\left(4-x^{2}\right)^{3}\cdot dx $$ $$=\int_{-2}^{2} {\left( 3x^{5}-24x^{3}-48^x-x^{6}+12x^{4}-48x^{2}+64 \right) dx}$$ (6.4.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$= \left( \frac{3x^{6}}{6}-\frac{24x^{4}}{4}-\frac{48x^{2}}{2}-\frac{x^{7}}{7}+\frac{12x^{5}}{5}-\frac{48x^{3}}{3}+64x \right)_{-2}^{2}=\color{blue}\frac{4096}{35} $$     (6.4.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \overrightarrow{q}=q_{x}\overrightarrow{i}+q_{y}\overrightarrow{j} $$     (6.4.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \overrightarrow{q}=\left({3x^{2}y+y^3}\right) \overrightarrow{i}+\left({3x+y^3}\right)\overrightarrow{j} $$     (6.4.13)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$The\; normal\; to\; the\; curved\; part\; of\; domain\; is\; \overrightarrow{j}\; and\; the\; normal\; to\; the\; straight\; line\; part\; is\; \overrightarrow{-j}$$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \oint{\overrightarrow{q}\cdot \overrightarrow{n}d\Gamma}=\int_{-2}^{2}- \left({3x^{2}y+y^3}\right) dx+ \int_{-2}^{2} \left({3x^{2}y+y^3}\right) dx $$ (6.4.14)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$For\; the\; straight\; line\; y=0\; and\; for\; curved\; part\; it\; is\; y=4-x^2$$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \oint{\overrightarrow{q}\cdot \overrightarrow{n}d\Gamma}=0+ \int_{-2}^{2} 3x^{2} {\left( 4-x^2 \right) }  +  \left(4-x^{2} \right)^{3}  dx = \color{blue}\frac{4096}{35} $$     (6.4.15)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$\color{blue}we\;have\;verified\;the\;divergence\;theorem$$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Text book (Fish and Belytschko) Problem 6.4.3

 * {| style="width:100%" border="0"

$$ Using\; the\; divergence\; theorem\; prove\; \color{blue}{\mathbf{ \oint_{\Gamma} nd\Gamma = 0 }} $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Solution

 * {| style="width:100%" border="0"

$$ We\;know\;that\; \color{blue}\oint_{\Gamma} \overrightarrow{q} \cdot \overrightarrow{n} d\Gamma= \int_{\Omega}\overrightarrow{\nabla}\cdot \overrightarrow{q}d\Omega $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ Assume\; \overrightarrow{q}\; is\; an\; unit\; vector,\; which\; means\; the\; vector\; fields\; are\; q_x=1\;,q_y=1\;. $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \overrightarrow{q} \cdot \overrightarrow{n}={n} $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \oint_{\Gamma} \overrightarrow{q} \cdot \overrightarrow{n} d\Gamma= \oint_{\Gamma} nd\Gamma $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \int_{\Omega}\overrightarrow{\nabla}\cdot \overrightarrow{q}d\Omega =\int_{\Omega}\left( \frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y} \right) d\Omega=\int_{\Omega}\left( \frac{\partial 1}{\partial x}+\frac{\partial 1}{\partial y} \right) d\Omega= \int_{\Omega}0\cdot d\Omega=\color{blue}0 $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \color{blue}Hence\; proved $$
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Problem Description
Use 2D LIBF with n=2,4,6,8 to solve G2DM1.0/D1 for $${{u}^{h}}$$ until accuracy $${{10}^{-6}}$$ at center (x,y)=(0,0).

Given

For biunit square $$\Omega =\bar{\omega }=\square $$

PDE : $$div(K\centerdot grad(u))+f=\rho c\frac{\partial u}{\partial t}$$, where K=I(identity matrix), f=0, $$\frac{\partial u}{\partial t}=0$$

Essential Boundary Condition: g=2 on $$\partial \Omega $$

Nature Boundary Condition: none

Solution
(1)Find Basis Functions

For 1D case Basis function:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{L_{i,m}}(x) = \prod\limits_{\begin{array}{ccccccccccccccc} {k = 1} \\ {k \ne i} \end{array}}^m {\frac}

$$     (Eq 6.5.1)
 * <p style="text-align:right">
 * }

For 2D case:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{N_I}(x,y) = {L_{i,m}}(x) \cdot {L_{j,n}}(y),whereI = i + (j - 1)m

$$     (Eq 6.5.2)
 * <p style="text-align:right">
 * }

Here are two example for basis functions, N(x):

(A)n=4, I=1



(B)n=4, I=7



(2)Find Matrix equations

Considering:
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

m(w,u) + k(w,u) = f(w,u)

$$     (Eq 6.5.3)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

m(w,u) = \int\limits_\Omega {w\rho c\frac} d\Omega

$$     (Eq 6.5.4)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

k(w,u) = \int\limits_\Omega {\nabla wK\nabla u} d\Omega

$$     (Eq 6.5.5)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

f(w,u) = - \int_{\Gamma h} {whd{\Gamma _h}}  + \int\limits_\Omega  {wf} d\Omega

$$     (Eq 6.5.6)
 * <p style="text-align:right">
 * }

Because $$\frac = 0$$, $$f = 0 $$ We get
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

Kd = F

$$     (Eq 6.5.7) In Eq 6.5.7
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{K_{ij}} = \int_ {(\nabla N_I^e)I(\nabla N_J^e)d{w^e}}

$$     (Eq 6.5.8)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

F = 0

$$     (Eq 6.5.9)
 * <p style="text-align:right">
 * }

(3)Reconstruct $${K_f}$$ and $${F_f}$$

Considering the essential boundary condition is the outline of the biunit, one special treatment is needed. Reconstruct $${K_f}$$ and $${F_f}$$, so that we can solve for matrix d. Firstly, considering
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{L_{i,m}}({x_j}) = {\delta _{i,m}} = \left\{ {\begin{array}{ccccccccccccccc} {1,forj = i} \\ {0,forj \ne i} \end{array}} \right.

$$     (Eq 6.5.10 ) For all $${d_i}$$ belong to the nodes on the essential boundary condition, we can get the value directly:
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:95%" |
 * style="width:95%" |

{d_i} = 2,where(i \in ess.b.c)

$$     (Eq 6.5.11) For the unknown $${d_i}$$, which local in the middle area of the biunit, whose number is $${(n - 2)^2}$$, we need to reconstruct $${K_f}$$ and $${F_f}$$ to solve for the unknown d.
 * <p style="text-align:right">
 * }

Considering Eq. 6.5.7, we can get
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{K_f}{d_f} = {F_f}

$$     (Eq 6.5.12) Comments(1): Dimension of $${K_f}$$ is $${n^2} \times {(n - 2)^2}$$, dimension of $${d_f}$$ is $${(n - 2)^2} \times 1$$ , dimension of $${F_f}$$ is $${n^2} \times 1$$
 * <p style="text-align:right">
 * }

Comments(2):In Eq. 6.5.12, $${F_f}$$ has absorbed the essential boundary condition, it means all "2" in $${d}$$ has timed the corresponding $${K_{i,j}}$$, and these products make the $${-F_f}$$ Matrix. Take advantage of Eq. 6.5.12, we can get $${d_f}$$.

(4)Reunion matrix $$d$$

In this step, update the unknown $${d_i}$$ in the $${d}$$ with $${d_f}$$, so we get the complete matrix $${d}$$. So the trial solution is :
 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

{u^h} = \sum\limits_{i = 1}^ {{d_i} \times {N_i}(x,y)}

$$     (Eq 6.5.13)
 * <p style="text-align:right">
 * }

Then plot $${u^h}$$:

For three different cases: m=n=2, m=n=4, m=n=6, we get the same trial solution, $${u^h=2}$$, shown below:



Exact solution is $${u(x,y)=2}$$, so the error in center point(0,0) between $${u^h}$$ and $${u(x,y)}$$ is 0 for all the three cases.

Comments(3):

Although it's too complicated to find exact solution via math way, we still can find exact solution though physic analysis. Firstly, because $$\frac = 0$$, the temperature field is static. Secondly, because $$f=0$$, there is no heat exchange between this biunit and the environment. So the total amount of heat in this biunit does not change.

Besides， K matrix is identity matrix, so the biunit is isotropic material, so the temperature is equal everywhere in the biunit. Considering the essential boundary condition is $${g=2}$$, we can get the conclusion that, the exact solution is $${u(x,y)=2}$$

Problem Description
Solve the general 2-D model with data set 1 (G2DM1.0/D1) using 2-D Linear Lagrangian element basis function (LLEBF) by increasing the number of elements until an accuracy of $$\displaystyle 10^{-6}$$ is reached at the center $$\displaystyle (x,y)=(0,0)$$.

Integrate $$\displaystyle \mathbf{k^{e}}$$ in parent coordinates $$\left \{ \xi_{i} \right \}$$. Find $$\displaystyle \mu$$ to integrate $$\displaystyle \mathbf{k^{e}}$$ exactly with Gauss-Legendre quadrature (GLQ)

Given
For G2DM1.0/D1, the governing partial differential equation (PDE) and given conditions are,


 * $$\displaystyle \frac{\partial}{\partial x_{i}} \left[ \boldsymbol{\kappa_{ij}} \frac{\partial u}{\partial x_{j}} \right ] + f = \rho c \frac{\partial u}{\partial t} $$


 * where
 * The domain $$\displaystyle \Omega = \overline{\omega}$$ is a biunit square


 * $$\displaystyle \boldsymbol{\kappa_{ij}}=\boldsymbol{I}$$


 * $$\displaystyle f = 0$$


 * $$\displaystyle \bar{m} \frac{\partial u}{\partial t } = 0 $$ (static case)

So the PDE simplifies to:
 * $$\displaystyle \frac{\partial}{\partial x_{i}} \left[ \boldsymbol{I} \frac{\partial u}{\partial x_{j}} \right ] = 0$$

The natural boundary condition is not prescribed and the essential boundary condition is defined below:


 * $$\displaystyle g=2$$ on $$\displaystyle \partial u$$

Generating LLEBF
The 2-d Lagrangian interpolation function can be formed as follows
 * $$\displaystyle N_{I}^{e}(x,y)=L_{i,m}^{e}(x)\cdot L_{i,m}^{e}(y)$$

where
 * $$\displaystyle I=i+(j-1)m$$

For linear basis functions there are two nodes per element in along each dimension $$\displaystyle n=m=2$$. So for each element there are four nodes and four corresponding basis functions.
 * $$\displaystyle N_{1}^{e}(x,y)=L_{1,2}^{e}(x)\cdot L_{1,2}^{e}(y)=\frac{(x-x_{2})(y-y_{2})}{(x_{1}-x_{2})(y_{1}-y_{2})}$$
 * $$\displaystyle N_{2}^{e}(x,y)=L_{2,2}^{e}(x)\cdot L_{1,2}^{e}(y)=\frac{(x-x_{1})(y-y_{2})}{(x_{2}-x_{1})(y_{1}-y_{2})}$$
 * $$\displaystyle N_{3}^{e}(x,y)=L_{1,2}^{e}(x)\cdot L_{2,2}^{e}(y)=\frac{(x-x_{2})(y-y_{1})}{(x_{1}-x_{2})(y_{2}-y_{1})}$$
 * $$\displaystyle N_{4}^{e}(x,y)=L_{2,2}^{e}(x)\cdot L_{2,2}^{e}(y)=\frac{(x-x_{1})(y-y_{1})}{(x_{2}-x_{1})(y_{2}-y_{1})}$$

For example if one element is used the basis function for the first element on that element is as follows:
 * $$\displaystyle N_{1}^{1}(x,y)=\frac{(x-1)(y-1)}{(-1-1)(-1-1)}$$

One can see from the plot of this basis function that it is equal to one at (-1,-1) and equal to zero at all other nodal points and thereby satisfies the CBS:



The second plot shows all four nodal basis functions for the single element plotted together. Again, each satisfies the CBS.

Solving using LLEBF
The following matlab code was adapted from Problem 6.5 to solve the (G2DM1.0/D1) using 2-D LLEBF.



Determining the Exact Solution
The exact solution is found simply by analyzing the data given in the problem. One can see the temperature field does not change with time as given by $$\displaystyle \frac = 0$$. There is no internal heat generation as given by $$\displaystyle f=0$$. Lastly the conductance matrix is equal to the identity matrix meaning that the material is isotropic. The temperature is given around the entire domain as 2 (essential boundary condition) and no flux is applied (natural boundary condition). So the biunit square is isotropic, has temperature two around its entire border with no flux, internal heat generation, or change in temperature.

The only possible conclusion is that the exact solution is there must be constant temperature across the entire domain equal to the temperature at the boundary: $$\displaystyle {u(x,y)=2}$$.

Alternatively we note that the PDE is in fact equal to the Laplace Equation. The solution $$\displaystyle {u(x,y)=2}$$ could also be obtained by solving the Laplace equation with the given boundary conditions.

So the temperature at the center is $$\displaystyle {u(0,0)=2}$$. The approximate solution yielded $$\displaystyle {u^h(0,0)=2}$$. We find the error to be zero using only one element.

Verifying the Solution Using CCX FEM Code
The solution was verified using Calculix with both uniform and distorted elements. In both cases the solution was again $$\displaystyle {u(x,y)=2}$$. The body was meshed using 8 node quadratic element.

code for uniform mesh
DEFINE CGX GEOMETRY<BR> DEFINE THE ELEMENTS THAT ENFORCE THE BOUNDARY ELEMENTS CCX INPUT FILE

CODE FOR NON UNIFORM MESH
CGX GEOMETRY CODE DEFINING THE NODES THAT ENFORCE THE BOUNDARY CONDITION CCX INPUT FILE





Integrating the Element Stiffness Matrix in Parent Coordinates with GLQ

 * $$\displaystyle \tilde{k}(w,u)=\int_{\Omega^{}}^{.} \nabla w \boldsymbol{\kappa}\nabla u d\Omega$$
 * $$\displaystyle \tilde{m}(w,u)=\int_{\Omega^{e}}^{.} w\rho c\frac{\partial u}{\partial t}d\Omega=0$$
 * $$\displaystyle \tilde{f}(w)=-\int_{\Gamma_{h}^{e}} whd\Gamma_{h}^{e}+\int_{\Omega}  wfd\Omega=0$$


 * $$\displaystyle \tilde{k}_{ij}^{e}=\int_{\omega^{e}}^{.} \nabla(N_{I}^{e}) \boldsymbol{\kappa}\nabla(N_{J}^{e})d\omega^{e}$$


 * $$\displaystyle \tilde{k}_{ij}^{e}=\int_{\omega^{e}}^{.} \left [\boldsymbol{J}^{e*-T}\nabla_\xi N_{I}^{e}(\xi))\right ] \boldsymbol{\kappa}(\xi) \left [\boldsymbol{J}^{e*-T}\nabla_\xi N_{J}^{e}(\xi))   \right ] |\boldsymbol{J}| d\omega^{e}$$


 * $$ \displaystyle {K^1}_{11} = \int_{-1}^1 \int_{-1}^1 \begin{bmatrix} (y-1)/4 &(x-1)/4 \end{bmatrix} \begin{bmatrix}1 &0 \\ 0 &1 \end{bmatrix} \begin{bmatrix} (y-1)/4 \\ (x-1)/4 \end{bmatrix} dx dy $$


 * $$ \displaystyle {K^1}_{11} = \int_{-1}^1 \int_{-1}^1 \frac{x^2}{16}+ \frac{x}{8}+ \frac{y^2}{16}+ \frac{y}{8} dx dy $$

The integrand is second order in both x and y. In order to integrate a second exactly with GLQ we must select at least $$\displaystyle \mu=2$$.


 * $$\displaystyle W_{1}=W_{2}=1, \xi_{1}= -\frac{1}{\sqrt{3} }, \xi_{2}=\frac{1}{\sqrt{3}}$$


 * $$\displaystyle I_{11}= 2 \left [W_{1}(\frac{x^2}{16}+ \frac{x}{8}+ \frac{y^2}{16}+ \frac{y}{8}) +W_{2}(\frac{x^2}{16}+ \frac{x}{8}+ \frac{y^2}{16}+ \frac{y}{8}) \right ] $$


 * $$I_{11}= 2 \left [W_{1}(\frac{\xi_1^2}{16}+ \frac{\xi_1}{8}+ \frac{\xi_1^2}{16}+ \frac{\xi_1}{8}) +W_{2}(\frac{\xi_2^2}{16}+ \frac{\xi_2}{8}+ \frac{\xi_2^2}{16}+ \frac{\xi_2}{8}) \right ] $$


 * $$\displaystyle I_{11}= 2/3$$

Similarly,
 * $$\displaystyle I_{12}= 1/3$$


 * $$\displaystyle I_{21}= 1/3$$


 * $$\displaystyle I_{22}= 1/3$$

Problem Description
Prove the following matrix identities


 * {| style="width:100%" border="0"

$$ \left( \mathbf{A \cdot B} \right)^T = \mathbf{B}^T \cdot \mathbf{A}^T $$     (6.7.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \left( \mathbf{A}^{-1} \right )^T = \left( \mathbf{A}^{T} \right )^{-1} = \mathbf{A}^{-T} $$     (6.7.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * for the cases when



\mathbf{A} = \begin{bmatrix} 1 & 1 & 1 \\ 2 & -1 & 3 \\ 3 & 2 & 6 \end{bmatrix} \qquad \text{and} \qquad \mathbf{B} = \begin{bmatrix} 1 & 3 & 5 \\ 1 & -4 & 1 \\ 2 & 5 & 8 \end{bmatrix} $$

Also explain the syntax of the WolframAlpha knowledge engine.

Solution

 * Solutions have been confirmed using WolframAlpha. See hyperlinks designated "WA".

Part 1: Prove equation 6.7.1
The transpose of a matrix involves writing the column values as rows and the row values as columns. In other words, $$\displaystyle [\mathbf{A}^T]_{ij} = [\mathbf{A}]_{ji}$$. Therefore,



\mathbf{A}^T = \begin{bmatrix} 1 & 2 & 3 \\ 1 & -1 & 2 \\ 1 & 3 & 6 \end{bmatrix} \qquad \text{and} \qquad \mathbf{B}^T = \begin{bmatrix} 1 & 1 & 2 \\ 3 & -4 & 5 \\ 5 & 1 & 8 \end{bmatrix} $$.

So,



\mathbf{B}^T \cdot \mathbf{A}^T = \begin{bmatrix} 1 & 1 & 2 \\ 3 & -4 & 5 \\ 5 & 1 & 8 \end{bmatrix} \cdot \begin{bmatrix} 1 & 2 & 3 \\ 1 & -1 & 2 \\ 1 & 3 & 6 \end{bmatrix} = \begin{bmatrix} (1+1+2) & (2-1+6) & (3+2+12) \\ (3-4+5) & (6+4+15) & (9-8+30) \\ (5+1+8) & (10-1+24) & (15+2+48) \end{bmatrix} = \begin{bmatrix} 4 & 7 & 17 \\ 4 & 25 & 31 \\ 14 & 33 & 65 \end{bmatrix} $$. [WA ]

Compare this with



(\mathbf{A \cdot B})^T = \left( \begin{bmatrix} 1 & 1 & 1 \\ 2 & -1 & 3 \\ 3 & 2 & 6 \end{bmatrix} \cdot \begin{bmatrix} 1 & 3 & 5 \\ 1 & -4 & 1 \\ 2 & 5 & 8 \end{bmatrix} \right)^T = \begin{bmatrix} (1+1+2) & (3-4+5) & (5+1+8) \\ (2-1+6) & (6+4+15) & (10-1+24) \\ (3+2+12) & (9-8+30) & (15+2+48) \end{bmatrix}^T = \begin{bmatrix} 4 & 7 & 17 \\ 4 & 25 & 31 \\ 14 & 33 & 65 \end{bmatrix} $$. [WA ]

Therefore,


 * {| style="width:10%" border="0"

$$\displaystyle \left( \mathbf{A \cdot B} \right)^T = \mathbf{B}^T \cdot \mathbf{A}^T $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }

Part 2: Prove equation 6.7.2
The first step in inverting a matrix is confirming that the determinant of the matrix does not equal zero. If $$\displaystyle \left | \mathbf{A} \right | = 0$$, then the matrix is not invertible.


 * For $$ \mathbf{A} = \begin{bmatrix}

1 & 1 & 1 \\ 2 & -1 & 3 \\ 3 & 2 & 6 \end{bmatrix} \text{,} \qquad \left | \mathbf{A} \right | = 1(-6-6) - 1(12-9) + 1(4+3) = -8 $$. [WA ]


 * For $$ \mathbf{B} = \begin{bmatrix}

1 & 3 & 5 \\ 1 & -4 & 1 \\ 2 & 5 & 8 \end{bmatrix} \text{,} \qquad \left | \mathbf{B} \right | = 1(-32-5) - 3(8-2) + 5(5+8) = 10 $$. [WA ]

Both matrices are invertible. Therefore,


 * $$ \mathbf{A}^{-1} = \begin{bmatrix}

1 & 1 & 1 \\ 2 & -1 & 3 \\ 3 & 2 & 6 \end{bmatrix}^{-1} = \frac{1}{8} \begin{bmatrix} 12 & 4 & -4 \\ 3 & -3 & 1 \\ -7 & -1 & 3 \end{bmatrix}$$[WA ]
 * $$\qquad \Rightarrow

\left( \mathbf{A}^{-1}\right)^T = \frac{1}{8} \begin{bmatrix} 12 & 3 & -7 \\ 4 & -3 & -1 \\ -4 & 1 & 3 \end{bmatrix}$$.


 * Similarly,
 * $$ \left( \mathbf{A}^{T}\right)^{-1} =

\begin{bmatrix} 1 & 2 & 3 \\ 1 & -1 & 2 \\ 1 & 3 & 6 \end{bmatrix}^{-1} = \frac{1}{8} \begin{bmatrix} 12 & 3 & -7 \\ 4 & -3 & -1 \\ -4 & 1 & 3 \end{bmatrix}$$. [WA ]


 * Therefore,


 * {| style="width:10%" border="0"

$$\displaystyle \left( \mathbf{A}^{-1} \right )^T = \left( \mathbf{A}^{T} \right )^{-1} = \mathbf{A}^{-T} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }

The same holds for matrix $$\displaystyle \mathbf{B}$$.


 * $$ \mathbf{B}^{-1} = \begin{bmatrix}

1 & 3 & 5 \\ 1 & -4 & 1 \\ 2 & 5 & 8 \end{bmatrix}^{-1} = \frac{1}{10} \begin{bmatrix} -37 & 1 & 23 \\ -6 & -2 & 4 \\ 13 & 1 & -7 \end{bmatrix}$$[WA ]
 * $$\qquad \Rightarrow

\left( \mathbf{B}^{-1}\right)^T = \frac{1}{10} \begin{bmatrix} -37 & -6 & 13 \\ 1 & -2 & 1 \\ 23 & 4 & -7 \end{bmatrix}$$.


 * Similarly,
 * $$ \left( \mathbf{B}^{T}\right)^{-1} =

\begin{bmatrix} 1 & 1 & 2 \\ 3 & -4 & 5 \\ 5 & 1 & 8 \end{bmatrix}^{-1} = \frac{1}{10} \begin{bmatrix} -37 & -6 & 13 \\ 1 & -2 & 1 \\ 23 & 4 & -7 \end{bmatrix}$$. [WA ]


 * Therefore,


 * {| style="width:10%" border="0"

$$\displaystyle \left( \mathbf{B}^{-1} \right )^T = \left( \mathbf{B}^{T} \right )^{-1} = \mathbf{B}^{-T} $$
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * style="width:10%; padding:10px; border:2px solid #8888aa" |
 * }

Both matrices confirm that the inverse and transpose matrix functions are commutative.

Part 3: Syntax of WolframAlpha
WolframAlpha utilizes the Mathematica core language. In general, the syntax follows the form $$ function[argument]$$. Note that the argument is contained within the square brackets. Let us consider two examples of linear algebra matrix functions - "transpose" and "inverse".
 * Also let us define the matrices $$\displaystyle \mathbf{M} = \begin{bmatrix}

1 & 2 \\ 6 & 12 \\ -60 & 2 \end{bmatrix} \qquad \text{and} \qquad \mathbf{N} = \begin{bmatrix} 1 & -2 \\ 6 & 12 \end{bmatrix} $$.


 * TRANSPOSE: To find the transpose of the matrices in WolframAlpha we enter the following:
 * For matrix $$\mathbf{M}$$ enter this: transpose[{{1,2},{6,12},{-60,2}}]
 * Each row is contained in single curly brackets. Each value in the row is separated by a comma. Each row is separated by a comma, and the entire matrix is contained in curly brackets.
 * WolframAlpha outputs: $$\mathbf{M}^T = \begin{bmatrix}

1 & 6 & -60 \\ 2 & 12 & 2 \end{bmatrix}$$. View webpage output here.


 * For matrix $$\mathbf{N}$$ enter this: transpose[{{1,-2},{6,12}}]
 * WolframAlpha outputs: $$\mathbf{N}^T = \begin{bmatrix}

1 & 6 \\ -2 & 12 \end{bmatrix}$$. View webpage output here.


 * INVERSE: Only square matrices (the number of rows and columns are equal) with non-zero determinants are invertible. Therefore matrix $$\displaystyle \mathbf{N}$$ is invertible only. To find the inverse of a matrix we enter the following:


 * For matrix $$\mathbf{N}$$ enter this: inverse[{{1,-2},{6,12}}]
 * "inverse" can also be abbreviated as "inv" in WolframAlpha.
 * WolframAlpha outputs: $$\mathbf{N}^{-1} = \frac{1}{24} \begin{bmatrix}

12 & 2 \\ -6 & 1 \end{bmatrix}$$. View webpage output here.


 * COMBINATION: Now let us consider the inverse of a transposed matrix (see equation 6.7.2 above). To find the inverse of the transposed matrix $$\displaystyle \mathbf{N}$$, in other words $$ \left( \mathbf{N}^T \right)^{-1}$$, we enter the following:


 * Enter this: inverse[transpose[{{1,-2},{6,12}}]]
 * The "innermost" square bracket pair is calculated first (transpose). Then any remaining functions are subsequently determined as we work our way outward. Inverse is the next, and last, matrix function.
 * WolframAlpha outputs: $$ \left( \mathbf{N}^T \right)^{-1} = \frac{1}{24} \begin{bmatrix}

12 & -6 \\ 2 & 1 \end{bmatrix}$$. View the webpage output here.


 * ANOTHER EXAMPLE: The inverse of a transpose of an inverse of a transpose results in the original square matrix.
 * $$ \left\{ \left[ \left( \mathbf{N}^T \right)^{-1} \right]^T \right\}^{-1} = \mathbf{N}$$. See in WolframAlpha.

Part 1: Verification of Position of Gauss Points and Corresponding Weights
Verify table of $$ \displaystyle \left \{ \left (w_{i},x_{i} \right ),i=1,...,5 \right \} $$ against NIST Handbook and FB p.89

NIST Handbook: Reference: http://en.wikipedia.org/wiki/Gaussian_quadrature

FB p.89:

Solution
Using excel, the values from the NIST Handbook are calculated to be: Which is the same from table 4.1 from FB p.89

Part 2: Solution to Integrals with User Written Gauss Quadrature code and Checked with Analytical Solution
FB p.91 pbs 4.6: Use Gauss quadrature to obtain exact values for the following integrals. Verify by analytical integration a) ::{| style="width:100%" border="0" $$ \displaystyle \int_{0}^{4}\left ( x^{2}+1\right )dx $$       (6.8.1) b)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \int_{-1}^{1}\left ( \xi^{4}+2\xi^{2}\right )d\xi $$       (6.8.3) c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB CODE.
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

4.6: Use three-point Gauss quadrature to evaluate the following integrals. Compare to the analytical integral a) ::{| style="width:100%" border="0" $$ \displaystyle \int_{-1}^{1} \frac{\xi}{\xi^{2}+1}   d\xi $$        (6.8.3) b) ::{| style="width:100%" border="0" $$ \displaystyle \int_{-1}^{1} cos^{2}(\pi\zeta )  d\xi $$       (6.8.4) c) Write a MATLAB code that utilizes function gauss.m and performs Gauss integration. Check your manual calculations against the MATLAB CODE.  4.8: The integral $$ \int_{-1}^{1}  \left ( 2\xi^{3}+2 \right )  d\xi $$ can be integrated exactly using two-point Gauss quadrature. How is the accuracy affected if: a) one-point quadrature is employed; b) three-point quadrature is employed. Check your solution against MATLAB code.
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Solution
4.6
 * {| style="width:100%" border="0"

$$ \displaystyle 2n_{gp}-1=3\Rightarrow n_{gp}=2 $$       (6.8.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle \begin{bmatrix} 1 & 1\\ \xi_{1} & \xi_{2}\\ \xi^{2}_{1}& \xi^{2}_{2}\\ \xi^{3}_{1}& \xi^{3}_{2} \end{bmatrix} \begin{bmatrix} W_{1}\\ W_{2} \end{bmatrix}= \begin{bmatrix} 2\\ 0\\ \frac{2}{3}\\ 0 \end{bmatrix} $$       (6.8.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle W_{1}=W_{2}=1, \xi_{1}= -\frac{1}{\sqrt{3} }, \xi_{2}=\frac{1}{\sqrt{3}}, a=0, b=4, dx=2d\xi, l=4 $$       (6.8.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle x=\frac{1}{2}(a+b)+\frac{1}{2}\xi(b-a) $$       (6.8.9)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle x=\frac{1}{2}(0+4)+\frac{1}{2}\xi(4-0) $$       (6.8.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle x=2+2\xi $$       (6.8.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle I=\int_{-1}^{1}((2+2\xi)^{2}+1)2d\xi=2\int_{-1}^{1}(4\xi^{2}+8\xi+5)d\xi $$ $$ I= 2 \left [W_{1}(4\xi^{2} _{1} +8\xi_{1}+5) +W_{2}(4\xi^{2}_{2}+8\xi_{2}+5) \right ] $$ $$ I= 2 \left [(4 \left ( -\frac{1}{\sqrt{3} } \right ) ^{2} +8 \left (-\frac{1}{\sqrt{3}}  \right ) +5) + (4 \left ( \frac{1}{\sqrt{3} }  \right ) ^{2} +8 \left (\frac{1}{\sqrt{3}}  \right ) +5)  \right ] $$       (6.8.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

I=2*(\frac{8}{3}+10)=\frac{76}{3}=25.\overline{33}

$$     (6.8.13) This is the same answer using WolframAlpha  and has been verified in the MATLAB code b)
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle n_{gp}=3 $$       (6.8.14) The integral is evaluated from -1 to 1
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle W_{1}=\frac{8}{9}, W_{2}=W_{3}=\frac{5}{9}$$ $$ \xi_{1}=0, \xi_{2}=-\frac{\sqrt{15}}{5},\xi_{3}=\frac{\sqrt{15}}{5} $$       (6.8.15)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle I=\int_{-1}^{1}\left ( \xi^{4}+2\xi^{2} \right )d\xi $$ $$ I=W_{1}\left ( \xi_{1} ^{4}+2\xi_{1}^{2} \right )+W_{2}\left ( \xi_{2} ^{4}+2\xi_{2}^{2} \right )+W_{3}\left ( \xi_{3} ^{4}+2\xi_{3}^{2} \right ) $$ $$ I=0 + \frac{5}{9} \left ( \left (-\frac{\sqrt{15}}{5} \right ) ^{4}+2\left (-\frac{\sqrt{15}}{5}  \right )^{2} \right )+\frac{5}{9} \left ( \left (\frac{\sqrt{15}}{5}  \right ) ^{4}+2\left (\frac{\sqrt{15}}{5} \right )^{2} \right )$$ $$ I=\frac{5}{9} \left (2\frac{225}{625}  +4\frac{15}{25}  \right ) $$
 * style="width:10%; |
 * style="width:10%; |

(6.8.16)
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

I=\frac{26}{15}=1.7\overline{33}

$$     (6.8.17) This is the same answer using WolframAlpha and has been verified in the MATLAB code
 * <p style="text-align:right">
 * }

4.6 a)
 * {| style="width:100%" border="0"

$$ \displaystyle n_{gp}=3 $$       (6.8.18) The integral is evaluated from -1 to 1
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle W_{1}=\frac{8}{9}, W_{2}=W_{3}=\frac{5}{9}$$ $$ \xi_{1}=0, \xi_{2}=-\frac{\sqrt{15}}{5},\xi_{3}=\frac{\sqrt{15}}{5} $$       (6.8.19)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle I=\int_{-1}^{1}\left ( \frac{\xi}{\xi^{2}+1} \right )d\xi $$ $$ I=\frac{2}{2}\left (W_{1} \left ( \frac{\xi_{1}}{\xi_{1}^{2}+1} \right )+W_{2} \left ( \frac{\xi_{2}}{\xi_{2}^{2}+1} \right )+W_{3} \left ( \frac{\xi_{3}}{\xi_{3}^{2}+1} \right ) \right ) $$ $$ I=0+\frac{5}{9} \left ( \frac{\left ( -\frac{\sqrt{15}}{5} \right )}{\left ( -\frac{\sqrt{15}}{5}  \right )^{2}+1} \right )+\frac{5}{9} \left ( \frac{\left ( \frac{\sqrt{15}}{5}  \right )}{\left ( -\frac{\sqrt{15}}{5}  \right )^{2}+1} \right )  $$ (6.8.20)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

I=0

$$     (6.8.21) This is the same answer using WolframAlpha and has been verified in the MATLAB code
 * <p style="text-align:right">
 * }

4.8:

b)
 * {| style="width:100%" border="0"

$$ \displaystyle n_{gp}=3 $$       (6.8.22) The integral is evaluated from -1 to 1
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle W_{1}=\frac{8}{9}, W_{2}=W_{3}=\frac{5}{9}$$ $$ \xi_{1}=0, \xi_{2}=-\frac{\sqrt{15}}{5},\xi_{3}=\frac{\sqrt{15}}{5} $$       (6.8.23)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }
 * {| style="width:100%" border="0"

$$ \displaystyle I=\int_{-1}^{1}\left ( cos^{2}(\pi \zeta ) \right )d\xi $$ $$ I=\frac{2}{2}\left (W_{1} \left ( cos^{2}(\pi \xi_{1} ) \right )+W_{2} \left ( cos^{2}(\pi \xi_{2} ) \right )+W_{3} \left ( cos^{2}(\pi \xi_{3} ) \right ) \right ) $$ $$ I=\frac{8}{9}  \left (cos^{2}(0)  \right ) +  \frac{5}{9} cos^{2} \left (\pi \left (-\frac{\sqrt{15}}{5}   \right )  \right ) +\frac{5}{9} cos^{2} \left (\pi  \frac{\sqrt{15}}{5}  \right ) $$ (6.8.24)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$ \displaystyle
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

I= 1.52996164

$$     (6.8.25) The exact answer using WolframAlpha is I=1, both the approximate answer and the exact solution was verified using the MATLAB code
 * <p style="text-align:right">
 * }

4.9 Using the MATLAB code, the integral is calulated to be I=4 using the Gauss quadrature for one, two, and three point quadrature. There are no changes to the accuracy.

Problem Description
Consider the heat conduction problem depicted in Figure 6.9.1. The coordinates are given in meters. The conductivity is isotropic, with $$\mathbf{D}=k\left[ \begin{matrix} 1 & 0 \\   0 & 1  \\ \end{matrix} \right]$$, and $$k=5\text{W}{}^{\text{o}}{{\text{C}}^{\text{-1}}}$$. The temperature $$T=0$$ is prescribed along edges AB and AD. The heat fluxes $$\bar{q}=0$$ and $$\bar{q}=20\text{W}{{\text{m}}^{\text{-1}}}$$ are prescribed on edges BC and CD, respectively. A constant heat source $$s=6\text{W}{{\text{m}}^{\text{-2}}}$$ is applied over the plate. Consider this problem modeled with 16 quadrilateral finite elements as shown in Figure 6.9.2.

Solve this problem using the finite element code given in Section 12.5.



Solution
By using the matlab codes downloaded from download, The post-processing results for temperature and flux are shown in Figure 6.9.3 and Figure 6.9.4.

All the matlab codes for this problem can be found in the file folder 'ch08'.

The main program is given in heat2d.m file. It calls the following functions: include_flags, [K,f,d] = preprocessor, [ke, fe] = heat2Delem(e), [K,f] = assembly(K,f,e,ke,fe), f = src_and_flux(f), [d,f_E] = solvedr(K,f,d) and postprocessor(d).

In the function of [K,f,d] = preprocessor, input file, which defines material properties, mesh data, vector and matrices initializations, essential and natural conditions and mesh generation, is read.

In the function of [ke, fe] = heat2Delem(e), the element conductance matrix and nodal source vector are integrated using Gauss quadrature. And then assembled by the function of [K,f] = assembly(K,f,e,ke,fe).

The f = src_and_flux(f) function calculates the nodal boundary flux vector.

In the postprocessor(d) function, the results are plotted in Figure 6.9.3 and Figure 6.9.4.

Problem Description
In this example, we consider a manufactured solution of the form

$$T={{\left( r-a \right)}^{2}}={{x}^{2}}+{{y}^{2}}-2a\sqrt{{{x}^{2}}+{{y}^{2}}}+{{a}^{2}}\ ,$$

defined over the domain of a square plate with a hole as shown in the figure below. For heat equation with isotropic conductivity and k=1, the corresponding source term is given by

$$s=-{{\nabla }^{2}}T=2a\frac{1}{\sqrt{{{x}^{2}}+{{y}^{2}}}}-4\ .$$

The essential boundary conditions on $${{\Gamma }_{T}}$$ are

$$T\left( r=a \right)=0,\quad \quad T\left( x=\pm b,y \right)={{a}^{2}}+{{b}^{2}}+{{y}^{2}}-2a\sqrt{{{y}^{2}}+{{b}^{2}}}\ .$$

The natural boundary conditions on $${{\Gamma }_{q}}$$ are $$\left( \overline{q}=-k{{\mathbf{n}}^{T}}\nabla T \right)$$

$$\overline{q}\left( x,y=b \right)=-\frac{\partial T}{\partial y}\left( x,y=b \right)=2b\left( \frac{2a}{\sqrt{{{x}^{2}}+{{y}^{2}}}}-1 \right)\ ,$$

$$\overline{q}\left( x,y=-b \right)=-\frac{\partial T}{\partial y}\left( x,y=-b \right)=2b\left( 1-\frac{2a}{\sqrt{{{x}^{2}}+{{y}^{2}}}} \right)\ .$$

Refer to p205-207 of F&B for detailed problem description.

Solution
Assume $$a=0.25$$, $$b=1$$, C(-1,2/3), G(1,2/3).

Then the manufactured solution is:

$$T={{\left( r-0.25 \right)}^{2}}={{x}^{2}}+{{y}^{2}}-0.5\sqrt{{{x}^{2}}+{{y}^{2}}}+1/16$$

The heat source is given by:

$$s=-{{\nabla }^{2}}T=\frac{0.5}{\sqrt{{{x}^{2}}+{{y}^{2}}}}-4$$

The essential boundary condition are:

$$\begin{matrix} T(r=0.25)=0, & T(x=\pm 1,y)=\frac{17}{16}+{{y}^{2}}-\sqrt{{{y}^{2}}+1}. \\ \end{matrix}$$

The natural boundary condition are:

$$\bar{q}(x,y=1)=2\left( \frac{0.5}{\sqrt{{{x}^{2}}+{{y}^{2}}}}-1 \right)$$

$$\bar{q}(x,y=-1)=2\left( 1-\frac{0.5}{\sqrt{{{x}^{2}}+{{y}^{2}}}} \right)$$

Modify the matlab codes provided by F&B for example 8.3, so that they can satisfy the boundary conditions.

The 40-element mesh is shown in the following figure

The temperature distribution for the 40-element mesh is shown:

The heat flux for the 40-element mesh is shown:

The temperature along CG are shown below:

As we can see from the above figure, the error is much larger than that given in the textbook(F&B). This error might be due to the radius of the hole in the middle(we assume a=0.25, while in the textbook, the radius is around 0.3).

Problem Description
Consider a chimney constructed of two isotropic materials: dense concrete ($$k=2.0W{}^{\circ }{{C}^{-1}}$$) and bricks ($$k=0.9W{}^{\circ }{{C}^{-1}}$$). The temperature of the hot gases on the inside surface of the chimney is $$140{}^{\circ }C$$, whereas the outside is exposed to the surrounding air, which is at $$T=10{}^{\circ }C$$. The dimensions of the chimney (in meters) are shown below. For the analysis, exploit the symmetry and consider 1/8 of the chimney cross sectional area. Consider a mesh of eight elements as shown below. Determine the temperature and flux in the two materials. Analyze the problem with 2×2, 4×4 and 8×8 quadrilateral elements for 1/8 of the problem domain. A2×2 finite element mesh is shown in Figure below. Symmetry implies insulated boundary conditions on edges AD and BC. Note that element boundaries have to coincide with the interface between the concrete and bricks.

Solution
We can modify the matlab codes provided by F&B for example 8.3, so that they can satisfy the boundary conditions.

The mesh with 2×2, 4×4, 8×8 quadrilateral elements for 1/8 of the problem domain are shown below: The temperature and flux for 2×2 finite element mesh in the two elements are shown below:

solution_vector_d =

140.0000 140.0000  140.0000   10.0000   10.0000   10.0000   96.0445   87.8624   53.8161

reactions_vector =

26.0688  51.8123   70.3496  -54.5474  -82.3755  -11.3079

Heat flux at Gauss Points

Element 1:

Element 2:

Element 3:

Element 4

The temperature and flux for 4×4 finite element mesh in the two elements are shown below:

The temperature and flux for 8×8 finite element mesh in the two elements are shown below:

The function for 8×8 mesh

Jiang

Part 1: Develop Weak Form
Develop the Weak form of General 2-Dimensional Model:


 * {| style="width:100%" border="0"

$$\displaystyle div(\underline{K}\cdot gradu)+f=\rho c \frac{\partial u}{\partial t}$$ (6.10.1)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

with


 * {| style="width:100%" border="0"

$$\displaystyle \underline{K}=\underline{I},f(\underline{x},t)=0,\frac{\partial u}{\partial t}=0,g=2,h=3,H=.5,u_\infty =0.1$$ (6.10.2)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Part 1: Weak Form
G2DM2.0/D1 Gives the following:


 * {| style="width:100%" border="0"

$$\displaystyle div(\underline{I}\cdot gradu)=0$$ (6.10.3)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle \vec n(\Gamma_h)\cdot\underline{I}\cdot gradu|_{\Gamma_h}=h=3$$ (6.10.4)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle u(\Gamma_g)=g=2$$ (6.10.5)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle H[u(\Gamma_H)-u_\infty ]=-\vec n(\Gamma_H)\cdot\underline{I}\cdot gradu(\Gamma_H)$$ (6.10.6)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Where $$\Gamma_g$$ represents the location on the boundary of the essential boundary condition, $$\Gamma_h$$ represents the location of the natural boundary condition, and $$\Gamma_H$$ represents the location of Newton's Law of Cooling boundary condition. This is seen in the following diagram:

The first step to formulate the weak form is to multiply by a weighting function $$w(x,y)$$ and integrate over the domain of the problem, represented by $$\Omega$$:


 * {| style="width:100%" border="0"

$$\displaystyle \int w(\underline{x} ,\underline{y})div(\underline{I}\cdot gradu)d\Omega=0$$ (6.10.7)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Multiplying the natural and Newton's Law of Cooling Boundary conditions by the same weighting function:


 * {| style="width:100%" border="0"

$$\displaystyle w(\Gamma_h)[\underline{n}\cdot (\Gamma_h)\cdot \underline{I}\cdot gradu|_{\Gamma_h}-3]=0$$ (6.10.8)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle w(\Gamma_H)[H[u(\Gamma_H)-u_\infty]+ \underline{n} (\Gamma_h)\underline{I}\cdot gradu|_{\Gamma_H}]=0$$ (6.10.9)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Using integration by parts ($$u=w$$ and $$dv=div(\underline{I}\cdot gradu)d\Omega$$) as used here and the fundamental theorem of calculus gives:


 * {| style="width:100%" border="0"

$$\displaystyle \int \frac{\partial}{\partial x_i}(wI_{ij}\frac{\partial u}{\partial x_j})d\Omega-\int\frac{\partial w}{\partial x_i}I_{ij}\frac{\partial u}{\partial x_j}d\Omega=0$$ (6.10.10)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Note we have switched from vector notation to index notation/Einstein notation. Applying the divergence theorem on the first term:


 * {| style="width:100%" border="0"

$$\displaystyle \int \frac{\partial}{\partial x_i}(wI_{ij}\frac{\partial u}{\partial x_j})d\Omega=\int wn_iI_{ij}\frac{\partial }{\partial x_j}d(\partial\Omega)$$ (6.10.11)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Since $$\partial\Omega$$ is the boundary of the domain, it can be split up in


 * {| style="width:100%" border="0"

$$\displaystyle \int wn_iI_{ij}\frac{\partial }{\partial x_j}d(\partial\Omega)=\int_{\Gamma_g}wn_iI_{ij}\frac{\partial u}{\partial x_j}d\Gamma_g + \int_{\Gamma_h}wn_i I_{ij} \frac{\partial u}{\partial x_j}d\Gamma_h+\int_{\Gamma_H}w n_i I_{ij}\frac{\partial u}{\partial x_j}d\Gamma_H$$ (6.10.12)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Recall, to remove the unknown flux on the essential boundary condition the arbitrary weighting function will be set to zero on the essential boundary, $$w(\Gamma_g)=0$$. Additionally, plugging in the boundary equations with the weighting function gives the weak form:


 * {| style="width:100%" border="0"

$$ \displaystyle \int_{\Gamma_h}w(x,1)3d\Gamma_h+\int_{\Gamma_H}w(0,y).5[u(0,y)-.1]d\Gamma_H-\int_\Omega\frac{\partial w}{\partial x_i}I_{ij}\frac{\partial u}{\partial x_j}d\Omega=0
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |

$$     (Eq 6.10.13)
 * <p style="text-align:right">
 * }

Part 2: Semi Discrete Equations
To develop the semidiscrete version first form the following functions with the LIBF:


 * {| style="width:100%" border="0"

$$\displaystyle u\simeq u^h = \sum_{n=1}^{a*b}d_n\prod_{(j=1,j\neq n)}^a \prod_{(k=1,k\neq n)}^b\frac{(x-x_j)(y-y_k)}{(x_n-x_j)(y_n-y_k)}=\sum_{n=1}^a*d_n*N_n(x)$$ (6.10.14)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle w\simeq w^h = \sum_{m=1}^{a*b}c_m\prod_{(j=1,j\neq m)}^a \prod_{(k=1,k\neq m)}^b\frac{(x-x_j)(y-y_k)}{(x_m-x_j)(y_m-y_k)}=\sum_{m=1}^a*c_m*N_m(x)$$ (6.10.15)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Substituting these equations into the weak form:


 * {| style="width:100%" border="0"

$$\displaystyle \int_{\Gamma_h}\sum_{m}c_mN_m(x,1)3d\Gamma_h+\int_{\Gamma_H}\sum_{m}c_mN_m(0,y).5[\sum_{n}d_nN_n(0,y)-.1]d\Gamma_H-\int_\Omega\frac{\partial \sum_{m}c_mN_m}{\partial x_i}I_{ij}\frac{\partial \sum_{n}d_nN_n}{\partial x_j}d\Omega=0 $$ (6.10.16)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

Factoring out the constant $$c_m$$ and recognizing subscripts m,n,i, and j mean to sum over:


 * {| style="width:100%" border="0"

$$\displaystyle c_m\Bigg [\int_{\Gamma_h}N_m(x,1)3d\Gamma_h+\int_{\Gamma_H}N_m(0,y).5[d_nN_n(0,y)-.1]d\Gamma_H-\int_\Omega\frac{\partial N_m}{\partial x_i}I_{ij}\frac{\partial d_nN_n}{\partial x_j}d\Omega\Bigg]=0   $$ (6.10.17)
 * style="width:10%; |
 * style="width:10%; |
 * <p style="text-align:right">
 * }

We will now move to put this into matrix form similarly to this previous work. In order to do this the numbering of the nodes must begin on the essential boundary. This allows partitioning of the matrices as before, with essential and free components. By numbering the nodes on the essential boundary first, this means the basis functions, $$\displaystyle N_m,N_n$$ first functions will be on the essential boundary (m and n correspond to the node number). With this recognition we can form the following:


 * {| style="width:100%" border="0"

$$\displaystyle \vec C=\begin{bmatrix} \vec C_E \\ \vec C_F\end{bmatrix}=\begin{bmatrix} 0 \\ \vec C_F\end{bmatrix} $$     (6.10.18)
 * <p style="text-align:right">
 * }

The zero in equation comes from satisfying the constraint breaking solution where the coefficient at the essential boundary condition for the only nonzero basis function coefficient on the weighting function is zero. Hence the subscript E in equation 5.3.11 is for essential. Looking at the coefficients for uh:


 * {| style="width:100%" border="0"

$$\displaystyle \vec d=\begin{bmatrix} \vec d_E \\ \vec d_F\end{bmatrix}=\begin{bmatrix} g \\ \vec d_F\end{bmatrix} $$     (6.10.19)
 * <p style="text-align:right">
 * }

Likewise here, g is the value at the essential boundary condition. Since there is only one nonzero basis function here, the coefficient (since the basis function is one) must be equal to the essential boundary condition value. The subscript F can be thought of as free and will be determined.

Now defining the following


 * {| style="width:100%" border="0"

$$\displaystyle \tilde k(N_m,N_n)=\int_0^1\int_0^1(\frac{d}{dx}N_m)I_{ij}(\frac{d}{dx}N_n)dxdy-\int_{\Gamma_H}.5Nn(0,y)d\Gamma_H $$     (6.10.20)
 * <p style="text-align:right">
 * }

This shows how the far left integral term of equation 5.3.10 can be written in components of a matrix. Note the additional term that contains the constants $$\displaystyle d_n$$ which is what the K matrix will multiply. Additionally defining:


 * {| style="width:100%" border="0"

$$\displaystyle K_{EE}=[\tilde k(N_{m=1\rightarrow a+b-1},N_{n=1\rightarrow a+b-1})] $$     (6.10.21)
 * <p style="text-align:right">
 * }

Where m increments from 1 to a+b-1 (a is number of x nodes and b is number of y nodes) and n increments from 1 to a+b-1. These end values (a+b-1) correspond to the number of nodes that are found on the essential boundary, hence being part of the $$\displaystyle K_{EE}$$ matrix. Continuing on with the defining:


 * {| style="width:100%" border="0"

$$\displaystyle K_{EF}=[\tilde k(N_{m=1\to a+b-1},N_{n=a+b\to a*b}] $$     (6.10.22)
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle K_{FE}=[\tilde k(N_{m=a+b\to a*b},N_{n=1\to a+b-1})] $$     (6.10.23)
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle K_{FF}=[\tilde k(N_{m=a+b\to a*b},N_{n=a+b\to a*b})] $$     (6.10.24)
 * <p style="text-align:right">
 * }

Giving:


 * {| style="width:100%" border="0"

$$\displaystyle \mathbf K=\begin{bmatrix} \mathbf K_{EE} \mathbf K_{EF} \\ \mathbf K_{FE} \mathbf K_{FF}\end{bmatrix} $$     (6.10.25)
 * <p style="text-align:right">
 * }

It can be seen that the subscripts E correspond the the essential nodes, and consequently the basis functions on the essential nodes (1 to a+b-1 where minus one takes into account double counting of the point shared with the connecting essential boundary borders). The F subscript corresponds to the free nodes (a+b to a*b) or free basis functions where the coefficient must still be determined.

Lastly, combining the constants into:
 * {| style="width:100%" border="0"

$$\displaystyle \mathbf F = [\tilde f_i] $$     (6.10.26)
 * <p style="text-align:right">
 * }

Where,


 * {| style="width:100%" border="0"

$$\displaystyle \tilde f_i=\int_{\Gamma_h}3N_m(x,1)d\Gamma_h - \int_{\Gamma_H}.5N_n(0,y)d\Gamma_H $$     (6.10.27)
 * <p style="text-align:right">
 * }

Putting it all together:


 * {| style="width:100%" border="0"

$$\displaystyle \vec C \cdot \mathbf K\mathbf d=\begin{bmatrix}0\\\vec C_F\end{bmatrix}\cdot \Big (\begin{bmatrix} \mathbf K_{EE} \mathbf K_{EF} \\ \mathbf K_{FE} \mathbf K_{FF}\end{bmatrix}\begin{bmatrix}\mathbf g\\\vec d_F\end{bmatrix}-\begin{bmatrix}\mathbf F_E\\ \mathbf F_F\end{bmatrix} \Big )=0 $$     (6.10.28) Now, our weighting function is arbitrary so what is in parenthesis corresponding to CF must equal zero. Writing this out explicitly as a simple example for four Nodes gives the following:
 * <p style="text-align:right">
 * }


 * {| style="width:100%" border="0"

$$\displaystyle \begin{bmatrix}K_{11}K_{12}K_{13}K_{14}\\K_{21}K_{22}K_{23}K_{24}\\K_{31}K_{32}K_{33}K_{34}\\K_{41}K_{42}K_{43}K_{44}\end{bmatrix}\begin{bmatrix}g\\g\\g\\d_4\end{bmatrix}=\begin{bmatrix}f_1+r_1\\f_2+r_2\\f_3+r_3\\f_4\end{bmatrix} $$     (6.10.29)
 * <p style="text-align:right">
 * }

The first 3 g's correspond to the essential boundary, so in order to solve for our unknown d4, in terms of our knowns, we rearrange to get the following:


 * {| style="width:100%" border="0"

$$\displaystyle \begin{bmatrix}K_{44}\end{bmatrix}\begin{bmatrix}d_4\end{bmatrix}=\begin{bmatrix}f_4-K_{41}g-K_{42}g-K_{43}g\end{bmatrix} $$     (6.10.30)
 * <p style="text-align:right">
 * }

Part 3: MATLAB Solution and Plot
The following is the MATLAB code used to find the coefficients using Part 2's steps:

This gave the following plot with 2 nodes in the x and y direction (4 nodes total) and a Temperature at x=y=.5 of 2.6437:



For 3 nodes in the x and y direction (9 nodes total) and a Temperature at x=y=.5 of 4.7215:



Lastly, 4 nodes in the x and y direction (16 nodes total) and a Temperature at x=y=.5 of 7.0725:



Comments
We were not able to become familiar enough with a FEA package to examine the accuracy of the solution, however examination of the plots can still be done. First the essential boundary appears to be held, but not on the correct boundaries. We were not able to correct this and believe the error is in the syntax of plotting, for examination of the output for the temperature shows that the essential boundary condition is held in the correct locations:

MATLAB OUTPUT FOR U^h

2.0000   3.3597    4.4995    5.4197    6.1200    6.6006    6.8614    6.9024    6.7237    2.0000    3.3456    4.4742    5.3856    6.0799    6.5571    6.8172    6.8602    6.6861    2.0000    3.3259    4.4383    5.3371    6.0223    6.4940    6.7521    6.7966    6.6276    2.0000    3.3006    4.3920    5.2742    5.9473    6.4112    6.6660    6.7117    6.5481    2.0000    3.2695    4.3351    5.1969    5.8548    6.3089    6.5590    6.6054    6.4478    2.0000    3.2328    4.2678    5.1052    5.7449    6.1868    6.4311    6.4777    6.3266    2.0000    3.1903    4.1900    4.9991    5.6174    6.0452    6.2823    6.3287    6.1844    2.0000    3.1422    4.1017    4.8785    5.4726    5.8839    6.1125    6.1583    6.0214    2.0000    3.0885    4.0030    4.7436    5.3102    5.7029    5.9217    5.9666    5.8375    2.0000    3.0290    3.8937    4.5942    5.1304    5.5024    5.7101    5.7535    5.6327    2.0000    2.9639    3.7740    4.4304    4.9332    5.2822    5.4775    5.5190    5.4069    2.0000    2.8930    3.6438    4.2522    4.7184    5.0423    5.2239    5.2632    5.1603    2.0000    2.8165    3.5031    4.0596    4.4862    4.7828    4.9495    4.9861    4.8928    2.0000    2.7343    3.3519    3.8526    4.2366    4.5037    4.6540    4.6876    4.6043    2.0000    2.6465    3.1902    3.6312    3.9694    4.2049    4.3377    4.3677    4.2950    2.0000    2.5529    3.0181    3.3954    3.6849    3.8866    4.0004    4.0265    3.9647    2.0000    2.4537    2.8354    3.1451    3.3828    3.5485    3.6422    3.6639    3.6136    2.0000    2.3488    2.6423    2.8805    3.0633    3.1908    3.2631    3.2800    3.2415    2.0000    2.2382    2.4387    2.6014    2.7263    2.8135    2.8630    2.8747    2.8486    2.0000    2.1219    2.2246    2.3079    2.3719    2.4166    2.4420    2.4480    2.4347    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000

Columns 10 through 18

6.3252   5.7069    4.8689    3.8111    2.5335    1.0362   -0.6810   -2.6178   -4.7745    6.2948    5.6865    4.8611    3.8185    2.5589    1.0821   -0.6118   -2.5228   -4.6508    6.2449    5.6487    4.8390    3.8156    2.5787    1.1283   -0.5358   -2.4134   -4.5046    6.1754    5.5936    4.8026    3.8025    2.5932    1.1747   -0.4529   -2.2897   -4.3356    6.0864    5.5211    4.7520    3.7790    2.6021    1.2213   -0.3633   -2.1517   -4.1441    5.9778    5.4313    4.6870    3.7451    2.6055    1.2682   -0.2668   -1.9995   -3.9299    5.8496    5.3240    4.6078    3.7010    2.6035    1.3154   -0.1634   -1.8329   -3.6930    5.7018    5.1994    4.5144    3.6465    2.5960    1.3627   -0.0533   -1.6520   -3.4335    5.5345    5.0575    4.4066    3.5818    2.5830    1.4103    0.0637   -1.4569   -3.1514    5.3476    4.8982    4.2846    3.5067    2.5646    1.4582    0.1875   -1.2474   -2.8466    5.1411    4.7215    4.1483    3.4213    2.5406    1.5062    0.3181   -1.0237   -2.5192    4.9150    4.5275    3.9977    3.3256    2.5112    1.5545    0.4556   -0.7857   -2.1692    4.6694    4.3161    3.8328    3.2196    2.4763    1.6031    0.5999   -0.5333   -1.7965    4.4042    4.0874    3.6537    3.1032    2.4360    1.6519    0.7510   -0.2667   -1.4012    4.1195    3.8413    3.4603    2.9766    2.3901    1.7009    0.9089    0.0142   -0.9832    3.8152    3.5778    3.2526    2.8396    2.3388    1.7502    1.0737    0.3095   -0.5426    3.4913    3.2970    3.0306    2.6923    2.2820    1.7996    1.2453    0.6190   -0.0794    3.1478    2.9988    2.7944    2.5347    2.2197    1.8494    1.4237    0.9428    0.4065    2.7848    2.6832    2.5439    2.3668    2.1519    1.8994    1.6090    1.2809    0.9150    2.4022    2.3503    2.2791    2.1885    2.0787    1.9496    1.8011    1.6333    1.4462    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000    2.0000

Columns 19 through 21

-7.1509  -9.7471  -12.5631   -6.9960   -9.5584  -12.3378   -6.8093   -9.3276  -12.0595   -6.5907   -9.0550  -11.7283   -6.3403   -8.7403  -11.3442   -6.0579   -8.3837  -10.9072   -5.7438   -7.9852  -10.4172   -5.3977   -7.5447   -9.8743   -5.0198   -7.0622   -9.2785   -4.6101   -6.5378   -8.6298   -4.1685   -5.9714   -7.9281   -3.6950   -5.3631   -7.1735   -3.1896   -4.7128   -6.3659   -2.6525   -4.0205   -5.5054   -2.0834   -3.2863   -4.5920   -1.4825   -2.5102   -3.6257   -0.8497   -1.6921   -2.6064   -0.1851   -0.8320   -1.5342    0.5114    0.0701   -0.4091    1.2398    1.0141    0.7690    2.0000    2.0000    2.0000

Additionally analyzing the natural and convective (Newton's law of cooling) boundary conditions in comparison to the essential boundary condition physically makes sense. Looking at the plot:



On the end where there shows a significant drop in temperature is the natural boundary condition. The natural boundary condition designates a flux from the material. This flux gives $$\displaystyle \frac{\partial u}{\partial x_i} = 3$$ which is the driving force of the energy loss. Note, this is as fast as the material can conduct the heat lost to the surroundings. On the convective boundary, the heat loss is slowed. The temperature difference is the driving condition on this boundary as well, and is similar in size $$u(0,y)-.1$$ where u at x=0 could be considered close to 2, the essential boundary value. However, there is the convective coefficient .5 multiplying this temperature difference, thus slowing the heat loss. This could allow for an insulation on this border, causing the increase in temperature seen in the plot. Also the larger temperature difference.