<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="221" HEIGHT="56" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="223" HEIGHT="57" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img1.png"
ALT="$\displaystyle \rho \frac{\partial^2 \vec u}{\partial t^2} + c \frac{\partial \vec u}{\partial t} - \div ( C \varepsilon(\vec u)) = \vec f$"> in <IMG
- WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img2.png"
- ALT="$ \Omega$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Omega$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
$\vec u=\vec u (\vec x,t)$
-->
<IMG
- WIDTH="79" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="81" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img4.png"
- ALT="$ \vec u=\vec u (\vec x,t)$"> is the deformation of the body, <IMG
- WIDTH="11" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \vec u=\vec u (\vec x,t)$">
+ is the deformation of the body, <IMG
+ WIDTH="13" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img5.png"
ALT="$ \rho$">
+
and <IMG
- WIDTH="10" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="11" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img6.png"
- ALT="$ c$"> the density and attenuation coefficient, and <IMG
- WIDTH="10" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ ALT="$ c$">
+ the density and attenuation coefficient, and <IMG
+ WIDTH="12" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img7.png"
- ALT="$ \vec f$"> external forces.
+ ALT="$ \vec f$">
+ external forces.
In addition, initial conditions
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="100" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="102" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img8.png"
ALT="$\displaystyle \vec u(\cdot, 0) = \vec u_0(\cdot)$"> on <IMG
- WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img2.png"
- ALT="$ \Omega$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Omega$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="48" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img9.png"
ALT="$\displaystyle \vec u(\vec x,t)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img10.png"
ALT="$\displaystyle = \vec d(\vec x,t) \qquad$"></TD>
<TD> </TD>
$\Gamma_D\subset\partial\Omega$
-->
<IMG
- WIDTH="67" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="68" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img11.png"
- ALT="$ \Gamma_D\subset\partial\Omega$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_D\subset\partial\Omega$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(3)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="96" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img12.png"
- ALT="$\displaystyle \vec n C \varepsilon(\vec u(\vec x,t))$"></TD>
+ ALT="$\displaystyle \vec n \ C \varepsilon(\vec u(\vec x,t))$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img13.png"
ALT="$\displaystyle = \vec b(\vec x,t) \qquad$"></TD>
<TD> </TD>
$\Gamma_N=\partial\Omega\backslash\Gamma_D$
-->
<IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img14.png"
- ALT="$ \Gamma_N=\partial\Omega\backslash\Gamma_D$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_N=\partial\Omega\backslash\Gamma_D$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img15.png"
ALT="$\displaystyle .$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
\vec u^T)$
-->
<IMG
- WIDTH="153" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="155" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img16.png"
ALT="$ \varepsilon(\vec u)= \tfrac 12 (\nabla \vec u + \nabla
-\vec u^T)$"> is the symmetric gradient of the displacement, also called the
+\vec u^T)$">
+ is the symmetric gradient of the displacement, also called the
<I>strain</I>. <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$"> is a tensor of rank 4, called the <I>stress-strain
+ ALT="$ C$">
+ is a tensor of rank 4, called the <I>stress-strain
tensor</I> that contains knowledge of the elastic strength of the material; its
symmetry properties make sure that it maps symmetric tensors of rank 2
(``matrices'' of dimension <IMG
- WIDTH="11" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="13" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img18.png"
- ALT="$ d$">, where <IMG
- WIDTH="11" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ ALT="$ d$">
+, where <IMG
+ WIDTH="13" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img18.png"
- ALT="$ d$"> is the spatial dimensionality) onto
+ ALT="$ d$">
+ is the spatial dimensionality) onto
symmetric tensors of the same rank. We will comment on the roles of the strain
and stress tensors more below. For the moment it suffices to say that we
interpret the term <!-- MATH
$\div ( C \varepsilon(\vec u))$
-->
<IMG
- WIDTH="85" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="87" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img19.png"
- ALT="$ \div ( C \varepsilon(\vec u))$"> as the vector with
+ ALT="$ \div ( C \varepsilon(\vec u))$">
+ as the vector with
components <!-- MATH
$\tfrac \partial{\partial x_j} C_{ijkl} \varepsilon(\vec u)_{kl}$
-->
<IMG
- WIDTH="103" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="105" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img20.png"
- ALT="$ \tfrac \partial{\partial x_j} C_{ijkl} \varepsilon(\vec u)_{kl}$">,
+ ALT="$ \tfrac \partial{\partial x_j} C_{ijkl} \varepsilon(\vec u)_{kl}$">
+,
where summation over indices <IMG
- WIDTH="38" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="39" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img21.png"
- ALT="$ j,k,l$"> is implied.
+ ALT="$ j,k,l$">
+ is implied.
<P>
The quasistatic limit of this equation is motivated as follows: each small
of the body. In general, this will be in the form of waves radiating away from
the location of the disturbance. Due to the presence of the damping term,
these waves will be attenuated on a time scale of, say, <IMG
- WIDTH="12" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img22.png"
- ALT="$ \tau$">. Now, assume
+ ALT="$ \tau$">
+. Now, assume
that all changes in external forcing happen on times scales that are
much larger than <IMG
- WIDTH="12" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img22.png"
- ALT="$ \tau$">. In that case, the dynamic nature of the change is
+ ALT="$ \tau$">
+. In that case, the dynamic nature of the change is
unimportant: we can consider the body to always be in static equilibrium,
i.e. we can assume that at all times the body satisfies
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="99" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img23.png"
ALT="$\displaystyle - \div ( C \varepsilon(\vec u))$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="27" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="29" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img24.png"
ALT="$\displaystyle = \vec f$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">in <IMG
- WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img2.png"
- ALT="$ \Omega$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Omega$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(5)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="48" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img9.png"
ALT="$\displaystyle \vec u(\vec x,t)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img10.png"
ALT="$\displaystyle = \vec d(\vec x,t) \qquad$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">on <IMG
- WIDTH="25" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="26" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img25.png"
- ALT="$ \Gamma_D$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_D$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(6)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="96" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img12.png"
- ALT="$\displaystyle \vec n C \varepsilon(\vec u(\vec x,t))$"></TD>
+ ALT="$\displaystyle \vec n \ C \varepsilon(\vec u(\vec x,t))$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img13.png"
ALT="$\displaystyle = \vec b(\vec x,t) \qquad$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">on <IMG
- WIDTH="25" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="26" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img26.png"
- ALT="$ \Gamma_N$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_N$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img15.png"
ALT="$\displaystyle .$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
$\vec f(\vec x,t)$
-->
<IMG
- WIDTH="45" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="47" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img27.png"
- ALT="$ \vec f(\vec x,t)$">.
+ ALT="$ \vec f(\vec x,t)$">
+.
<P>
While these equations are sufficient to describe small deformations, computing
large deformations is a little more complicated. To do so, let us first
introduce a stress variable <IMG
- WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img28.png"
- ALT="$ \sigma$">, and write the differential equations in
+ ALT="$ \sigma$">
+, and write the differential equations in
terms of the stress:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="52" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="54" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img29.png"
ALT="$\displaystyle - \div\sigma$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="27" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="29" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img24.png"
ALT="$\displaystyle = \vec f$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">in <IMG
- WIDTH="33" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="34" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img30.png"
- ALT="$ \Omega(t)$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Omega(t)$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(8)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="48" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img9.png"
ALT="$\displaystyle \vec u(\vec x,t)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img10.png"
ALT="$\displaystyle = \vec d(\vec x,t) \qquad$"></TD>
<TD> </TD>
$\Gamma_D\subset\partial\Omega(t)$
-->
<IMG
- WIDTH="85" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="86" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img31.png"
- ALT="$ \Gamma_D\subset\partial\Omega(t)$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_D\subset\partial\Omega(t)$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(9)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="96" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img12.png"
- ALT="$\displaystyle \vec n C \varepsilon(\vec u(\vec x,t))$"></TD>
+ ALT="$\displaystyle \vec n \ C \varepsilon(\vec u(\vec x,t))$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="98" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img13.png"
ALT="$\displaystyle = \vec b(\vec x,t) \qquad$"></TD>
<TD> </TD>
$\Gamma_N=\partial\Omega(t)\backslash\Gamma_D$
-->
<IMG
- WIDTH="115" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="116" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img32.png"
- ALT="$ \Gamma_N=\partial\Omega(t)\backslash\Gamma_D$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_N=\partial\Omega(t)\backslash\Gamma_D$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img15.png"
ALT="$\displaystyle .$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
Note that these equations are posed on a domain <IMG
- WIDTH="33" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="34" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img30.png"
- ALT="$ \Omega(t)$"> that
+ ALT="$ \Omega(t)$">
+ that
changes with time, with the boundary moving according to the
displacements <!-- MATH
$\vec u(\vec x,t)$
-->
<IMG
- WIDTH="48" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img33.png"
- ALT="$ \vec u(\vec x,t)$"> of the points on the boundary. To
+ ALT="$ \vec u(\vec x,t)$">
+ of the points on the boundary. To
complete this system, we have to specify the relationship between the
stress and the strain, as follows:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="81" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="82" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img34.png"
ALT="$\displaystyle \dot\sigma = C \varepsilon (\dot{\vec u}),$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
where a dot indicates a time derivative. Both the stress <IMG
- WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img28.png"
- ALT="$ \sigma$"> and the
+ ALT="$ \sigma$">
+ and the
strain <!-- MATH
$\varepsilon(\vec u)$
-->
<IMG
- WIDTH="33" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="34" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img35.png"
- ALT="$ \varepsilon(\vec u)$"> are symmetric tensors of rank 2.
+ ALT="$ \varepsilon(\vec u)$">
+ are symmetric tensors of rank 2.
<P>
Numerically, this system is solved as follows: first, we discretize
the time component using a backward Euler scheme. This leads to a
discrete equilibrium of force at time step <IMG
- WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img36.png"
- ALT="$ n$">:
+ ALT="$ n$">
+:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="105" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="106" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img37.png"
ALT="$\displaystyle -\div\sigma^n = f^n,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="165" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="167" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img38.png"
ALT="$\displaystyle \sigma^n = \sigma^{n-1} + C \varepsilon (\Delta \vec u^n),$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
$\Delta \vec u^n$
-->
<IMG
- WIDTH="35" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img39.png"
- ALT="$ \Delta \vec u^n$"> the incremental displacement for time step
+ ALT="$ \Delta \vec u^n$">
+ the incremental displacement for time step
<IMG
- WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img36.png"
- ALT="$ n$">. This way, if we want to solve for the displacement increment, we
+ ALT="$ n$">
+. This way, if we want to solve for the displacement increment, we
have to solve the following system:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="107" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="108" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img40.png"
ALT="$\displaystyle - \div C \varepsilon(\Delta\vec u^n)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="108" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="110" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img41.png"
ALT="$\displaystyle = \vec f + \div\sigma^{n-1}$"></TD>
<TD> </TD>
$\Omega(t_{n-1})$
-->
<IMG
- WIDTH="58" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="59" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img42.png"
- ALT="$ \Omega(t_{n-1})$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Omega(t_{n-1})$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(14)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="70" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="71" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img43.png"
ALT="$\displaystyle \Delta \vec u^n(\vec x,t)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="195" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="196" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img44.png"
ALT="$\displaystyle = \vec d(\vec x,t_n) - \vec d(\vec x,t_{n-1}) \qquad$"></TD>
<TD> </TD>
$\Gamma_D\subset\partial\Omega(t_{n-1})$
-->
<IMG
- WIDTH="110" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="111" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img45.png"
- ALT="$ \Gamma_D\subset\partial\Omega(t_{n-1})$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_D\subset\partial\Omega(t_{n-1})$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img3.png"
ALT="$\displaystyle ,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(15)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
- WIDTH="118" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="119" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img46.png"
- ALT="$\displaystyle \vec n C \varepsilon(\Delta \vec u^n(\vec x,t))$"></TD>
+ ALT="$\displaystyle \vec n \ C \varepsilon(\Delta \vec u^n(\vec x,t))$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
- WIDTH="195" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="196" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img47.png"
ALT="$\displaystyle = \vec b(\vec x,t_n)-\vec b(\vec x,t_{n-1}) \qquad$"></TD>
<TD> </TD>
$\Gamma_N=\partial\Omega(t_{n-1})\backslash\Gamma_D$
-->
<IMG
- WIDTH="140" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="141" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img48.png"
- ALT="$ \Gamma_N=\partial\Omega(t_{n-1})\backslash\Gamma_D$"><IMG
- WIDTH="7" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \Gamma_N=\partial\Omega(t_{n-1})\backslash\Gamma_D$">
+<IMG
+ WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img15.png"
ALT="$\displaystyle .$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
\{v\in H^1(\Omega(t_{n-1}))^d: v|_{\Gamma_D}=\vec d(\cdot,t_n) - \vec d(\cdot,t_{n-1})\}$
-->
<IMG
- WIDTH="393" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="394" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img49.png"
ALT="$ \Delta \vec u^n \in
\{v\in H^1(\Omega(t_{n-1}))^d: v\vert _{\Gamma_D}=\vec d(\cdot,t_n) - \vec d(\cdot,t_{n-1})\}$">
+
such that
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
$\vec b = 0$
-->
<IMG
- WIDTH="42" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="44" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img51.png"
- ALT="$ \vec b = 0$">, and that the deformation of the
+ ALT="$ \vec b = 0$">
+, and that the deformation of the
body is driven by body forces <IMG
- WIDTH="10" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="12" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img7.png"
- ALT="$ \vec f$"> and prescribed boundary displacements
+ ALT="$ \vec f$">
+ and prescribed boundary displacements
<IMG
- WIDTH="13" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img52.png"
- ALT="$ \vec d$"> alone. It is also worth noting that when integrating by parts, we
+ ALT="$ \vec d$">
+ alone. It is also worth noting that when integrating by parts, we
would get terms of the form <!-- MATH
$(C \varepsilon(\Delta\vec u^n), \nabla \varphi
)_{\Omega(t_{n-1})}$
-->
<IMG
- WIDTH="157" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="158" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img53.png"
ALT="$ (C \varepsilon(\Delta\vec u^n), \nabla \varphi
-)_{\Omega(t_{n-1})}$">, but that we replace it with the term involving the
+)_{\Omega(t_{n-1})}$">
+, but that we replace it with the term involving the
symmetric gradient <!-- MATH
$\varepsilon(\varphi)$
-->
<IMG
- WIDTH="33" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="35" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img54.png"
- ALT="$ \varepsilon(\varphi)$"> instead of <!-- MATH
+ ALT="$ \varepsilon(\varphi)$">
+ instead of <!-- MATH
$\nabla\varphi$
-->
<IMG
- WIDTH="27" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="28" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img55.png"
- ALT="$ \nabla\varphi$">. Due to
+ ALT="$ \nabla\varphi$">
+. Due to
the symmetry of <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$">, the two terms are equivalent, but the symmetric version
+ ALT="$ C$">
+, the two terms are equivalent, but the symmetric version
avoids a potential for round-off to render the resulting matrix slightly
non-symmetric.
<P>
The system at time step <IMG
- WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img36.png"
- ALT="$ n$">, to be solved on the old domain
+ ALT="$ n$">
+, to be solved on the old domain
<!-- MATH
$\Omega(t_{n-1})$
-->
<IMG
- WIDTH="58" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="59" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img42.png"
- ALT="$ \Omega(t_{n-1})$">, has exactly the form of a stationary elastic
+ ALT="$ \Omega(t_{n-1})$">
+, has exactly the form of a stationary elastic
problem, and is therefore similar to what we have already implemented
in previous example programs. We will therefore not comment on the
space discretization beyond saying that we again use lowest order
$\sigma^{n-1}$
-->
<IMG
- WIDTH="38" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="39" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img56.png"
- ALT="$ \sigma^{n-1}$"> to compute the next incremental
+ ALT="$ \sigma^{n-1}$">
+ to compute the next incremental
displacement, i.e. we need to compute it at the end of the time step
to make sure it is available for the next time step. Essentially,
the stress variable is our window to the history of deformation of
<P>
As indicated above, we need to have the stress variable <IMG
- WIDTH="21" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="23" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img57.png"
- ALT="$ \sigma^n$"> available
+ ALT="$ \sigma^n$">
+ available
when computing time step <IMG
- WIDTH="40" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="41" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img58.png"
- ALT="$ n+1$">, and we can compute it using
+ ALT="$ n+1$">
+, and we can compute it using
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="165" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="167" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img59.png"
ALT="$\displaystyle \sigma^n = \sigma^{n-1} + C \varepsilon (\Delta \vec u^n).$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<BR CLEAR="ALL"><P></P>
There are, despite the apparent simplicity of this equation, two questions
that we need to discuss. The first concerns the way we store <IMG
- WIDTH="21" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="23" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img57.png"
- ALT="$ \sigma^n$">: even
+ ALT="$ \sigma^n$">
+: even
if we compute the incremental updates <!-- MATH
$\Delta\vec u^n$
-->
<IMG
- WIDTH="35" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img39.png"
- ALT="$ \Delta \vec u^n$"> using lowest-order
+ ALT="$ \Delta \vec u^n$">
+ using lowest-order
finite elements, then its symmetric gradient <!-- MATH
$\varepsilon(\Delta\vec u^n)$
-->
<IMG
- WIDTH="55" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="56" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img60.png"
- ALT="$ \varepsilon(\Delta\vec u^n)$"> is
+ ALT="$ \varepsilon(\Delta\vec u^n)$">
+ is
in general still a function that is not easy to describe. In particular, it is
not a piecewise constant function, and on general meshes (with cells that are
not rectangles parallel to the coordinate axes) or with non-constant
stress-strain tensors <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$"> it is not even a bi- or trilinear function. Thus, it
+ ALT="$ C$">
+ it is not even a bi- or trilinear function. Thus, it
is a priori not clear how to store <IMG
- WIDTH="21" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="23" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img57.png"
- ALT="$ \sigma^n$"> in a computer program.
+ ALT="$ \sigma^n$">
+ in a computer program.
<P>
To decide this, we have to see where it is used. The only place where we
$(\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}$
-->
<IMG
- WIDTH="134" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="135" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img61.png"
- ALT="$ (\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}$">. In practice, we of
+ ALT="$ (\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}$">
+. In practice, we of
course replace this term by numerical quadrature:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="523" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="518" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img62.png"
ALT="$\displaystyle (\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})} = \sum_{K\s...
-...thbb{T}}} \sum_q w_q \sigma^{n-1}(\vec x_q) : \varepsilon(\varphi(\vec x_q)),$"></TD>
+...athbb{T}}} \sum_q w_q \ \sigma^{n-1}(\vec x_q) : \varepsilon(\varphi(\vec x_q),$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(19)</TD></TR>
</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
where <IMG
- WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="23" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img63.png"
- ALT="$ w_q$"> are the quadrature weights and <IMG
- WIDTH="20" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ w_q$">
+ are the quadrature weights and <IMG
+ WIDTH="21" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img64.png"
- ALT="$ \vec x_q$"> the quadrature points on
+ ALT="$ \vec x_q$">
+ the quadrature points on
cell <IMG
- WIDTH="18" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img65.png"
- ALT="$ K$">. This should make clear that what we really need is not the stress
+ ALT="$ K$">
+. This should make clear that what we really need is not the stress
<!-- MATH
$\sigma^{n-1}$
-->
<IMG
- WIDTH="38" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="39" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img56.png"
- ALT="$ \sigma^{n-1}$"> in itself, but only the values of the stress in the quadrature
+ ALT="$ \sigma^{n-1}$">
+ in itself, but only the values of the stress in the quadrature
points on all cells. This, however, is a simpler task: we only have to provide
a data structure that is able to hold one symmetric tensor of rank 2 for each
quadrature point on all cells (or, since we compute in parallel, all
$\varepsilon(\Delta \vec u^n(\vec x_q))$
-->
<IMG
- WIDTH="84" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="85" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img66.png"
- ALT="$ \varepsilon(\Delta \vec u^n(\vec x_q))$">, multiply it by the stress-strain tensor <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ ALT="$ \varepsilon(\Delta \vec u^n(\vec x_q))$">
+, multiply it by the stress-strain tensor <IMG
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$">, and use the
+ ALT="$ C$">
+, and use the
result to update the stress <!-- MATH
$\sigma^n(\vec x_q)$
-->
<IMG
- WIDTH="50" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="52" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img67.png"
- ALT="$ \sigma^n(\vec x_q)$"> at quadrature point <IMG
- WIDTH="11" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \sigma^n(\vec x_q)$">
+ at quadrature point <IMG
+ WIDTH="12" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img68.png"
- ALT="$ q$">.
+ ALT="$ q$">
+.
<P>
The second complication is not visible in our notation as chosen above. It is
$\Delta u^n$
-->
<IMG
- WIDTH="34" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="36" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img69.png"
- ALT="$ \Delta u^n$"> on the domain <!-- MATH
+ ALT="$ \Delta u^n$">
+ on the domain <!-- MATH
$\Omega(t_{n-1})$
-->
<IMG
- WIDTH="58" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="59" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img42.png"
- ALT="$ \Omega(t_{n-1})$">,
+ ALT="$ \Omega(t_{n-1})$">
+,
and then use this displacement increment to both update the stress as well as
move the mesh nodes around to get to <!-- MATH
$\Omega(t_n)$
-->
<IMG
- WIDTH="41" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="43" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img70.png"
- ALT="$ \Omega(t_n)$"> on which the next increment
+ ALT="$ \Omega(t_n)$">
+ on which the next increment
is computed. What we have to make sure, in this context, is that moving the
mesh does not only involve moving around the nodes, but also making
corresponding changes to the stress variable: the updated stress is a variable
$\Delta\vec u$
-->
<IMG
- WIDTH="26" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="28" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img71.png"
- ALT="$ \Delta\vec u$"> can be decomposed into three parts, a linear translation (the constant part
+ ALT="$ \Delta\vec u$">
+ can be decomposed into three parts, a linear translation (the constant part
of the displacement field in the neighborhood of a point), a dilational
component (that part of the gradient of the displacement field that has a
nonzero divergence), and a rotation. A linear translation of the material does
$\Delta\vec u=(y, -x)^T$
-->
<IMG
- WIDTH="107" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="108" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img72.png"
- ALT="$ \Delta\vec u=(y, -x)^T$">, with which <!-- MATH
+ ALT="$ \Delta\vec u=(y, -x)^T$">
+, with which <!-- MATH
$\varepsilon(\Delta \vec u)=0$
-->
<IMG
- WIDTH="75" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="77" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img73.png"
- ALT="$ \varepsilon(\Delta \vec u)=0$">). Nevertheless, if the the material was pre-stressed in a certain
+ ALT="$ \varepsilon(\Delta \vec u)=0$">
+). Nevertheless, if the the material was pre-stressed in a certain
direction, then this direction will be rotated along with the material. To
this end, we have to define a rotation matrix <!-- MATH
$R(\Delta \vec u^n)$
-->
<IMG
- WIDTH="60" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="61" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img74.png"
- ALT="$ R(\Delta \vec u^n)$"> that
+ ALT="$ R(\Delta \vec u^n)$">
+ that
describes, in each point the rotation due to the displacement increments. It
is not hard to see that the actual dependence of <IMG
- WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img75.png"
- ALT="$ R$"> on <!-- MATH
+ ALT="$ R$">
+ on <!-- MATH
$\Delta \vec u^n$
-->
<IMG
- WIDTH="35" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img39.png"
- ALT="$ \Delta \vec u^n$"> can
+ ALT="$ \Delta \vec u^n$">
+ can
only be through the curl of the displacement, rather than the displacement
itself or its full gradient (as mentioned above, the constant components of
the increment describe translations, its divergence the dilational modes, and
the curl the rotational modes). Since the exact form of <IMG
- WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img75.png"
- ALT="$ R$"> is cumbersome, we
+ ALT="$ R$">
+ is cumbersome, we
only state it in the program code, and note that the correct updating formula
for the stress variable is then
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="298" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="299" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img76.png"
ALT="$\displaystyle \sigma^n = R(\Delta \vec u^n)^T [\sigma^{n-1} + C \varepsilon (\Delta \vec u^n)] R(\Delta \vec u^n).$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
A^K_{ij}$
-->
<IMG
- WIDTH="103" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="105" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img77.png"
ALT="$ A_{ij} = \sum_K
-A^K_{ij}$"> built up of local contributions on each cell <IMG
- WIDTH="18" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+A^K_{ij}$">
+ built up of local contributions on each cell <IMG
+ WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img65.png"
- ALT="$ K$"> with entries
+ ALT="$ K$">
+ with entries
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="170" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="171" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img78.png"
- ALT="$\displaystyle A^K_{ij} = (C \varepsilon(\varphi_i), \varepsilon(\varphi_j))_K;$"></TD>
+ ALT="$\displaystyle A^K_{ij} = (C \varepsilon(\varphi_j), \varepsilon(\varphi_i))_K;$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(21)</TD></TR>
</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
In practice, <IMG
- WIDTH="27" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="29" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img79.png"
- ALT="$ A^K$"> is computed using numerical quadrature according to the
+ ALT="$ A^K$">
+ is computed using numerical quadrature according to the
formula
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="275" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="277" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img80.png"
ALT="$\displaystyle A^K_{ij} = \sum_q w_q [\varepsilon(\varphi_i(\vec x_q)) : C : \varepsilon(\varphi_j(\vec x_q))],$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
with quadrature points <IMG
- WIDTH="20" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="21" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img64.png"
- ALT="$ \vec x_q$"> and weights <IMG
- WIDTH="21" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \vec x_q$">
+ and weights <IMG
+ WIDTH="23" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img63.png"
- ALT="$ w_q$">. We have built these
+ ALT="$ w_q$">
+. We have built these
contributions before, in step-8 and step-17, but in both of these cases we
have done so rather clumsily by using knowledge of how the rank-4 tensor <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
ALT="$ C$">
+
is composed, and considering individual elements of the strain tensors
<!-- MATH
$\varepsilon(\varphi_i),\varepsilon(\varphi_j)$
-->
<IMG
- WIDTH="83" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="84" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img81.png"
- ALT="$ \varepsilon(\varphi_i),\varepsilon(\varphi_j)$">. This is not really
+ ALT="$ \varepsilon(\varphi_i),\varepsilon(\varphi_j)$">
+. This is not really
convenient, in particular if we want to consider more complicated elasticity
models than the isotropic case for which <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$"> had the convenient form
+ ALT="$ C$">
+ had the convenient form
<!-- MATH
- $c_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl}
+ $C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl}
+ \delta_{il} \delta_{jk})$
-->
<IMG
- WIDTH="235" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="241" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img82.png"
- ALT="$ c_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl}
-+ \delta_{il} \delta_{jk})$">. While we in fact do not use a more complicated
+ ALT="$ C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl}
++ \delta_{il} \delta_{jk})$">
+. While we in fact do not use a more complicated
form than this in the present program, we nevertheless want to write it in a
way that would easily allow for this. It is then natural to introduce
classes that represent symmetric tensors of rank 2 (for the strains and
stresses) and 4 (for the stress-strain tensor <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$">). Fortunately, deal.II
+ ALT="$ C$">
+). Fortunately, deal.II
provides these: the <TT>SymmetricTensor<rank,dim></TT> class template
provides a full-fledged implementation of such tensors of rank <TT>rank</TT>
(which needs to be an even number) and dimension <TT>dim</TT>.
<P>
What we then need is two things: a way to create the stress-strain rank-4
tensor <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$"> as well as to create a symmetric tensor of rank 2 (the strain
+ ALT="$ C$">
+ as well as to create a symmetric tensor of rank 2 (the strain
tensor) from the gradients of a shape function <IMG
- WIDTH="19" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="20" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img83.png"
- ALT="$ \varphi_i$"> at a quadrature
+ ALT="$ \varphi_i$">
+ at a quadrature
point <IMG
- WIDTH="20" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="21" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img64.png"
- ALT="$ \vec x_q$"> on a given cell. At the top of the implementation of this
+ ALT="$ \vec x_q$">
+ on a given cell. At the top of the implementation of this
example program, you will find such functions. The first one,
<TT>get_stress_strain_tensor</TT>, takes two arguments corresponding to
the Lamé constants <IMG
- WIDTH="12" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img84.png"
- ALT="$ \lambda$"> and <IMG
- WIDTH="13" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ \lambda$">
+ and <IMG
+ WIDTH="14" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img85.png"
- ALT="$ \mu$"> and returns the stress-strain tensor
+ ALT="$ \mu$">
+ and returns the stress-strain tensor
for the isotropic case corresponding to these constants (in the program, we
will choose constants corresponding to steel); it would be simple to replace
this function by one that computes this tensor for the anisotropic case, or
taking into account crystal symmetries, for example. The second one,
<TT>get_strain</TT> takes an object of type <TT>FEValues</TT> and indices
<IMG
- WIDTH="9" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="10" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img86.png"
- ALT="$ i$"> and <IMG
- WIDTH="11" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ ALT="$ i$">
+ and <IMG
+ WIDTH="12" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img68.png"
- ALT="$ q$"> and returns the symmetric gradient, i.e. the strain,
+ ALT="$ q$">
+ and returns the symmetric gradient, i.e. the strain,
corresponding to shape function <!-- MATH
$\varphi_i(\vec x_q)$
-->
<IMG
- WIDTH="48" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img87.png"
- ALT="$ \varphi_i(\vec x_q)$">, evaluated on the cell
+ ALT="$ \varphi_i(\vec x_q)$">
+, evaluated on the cell
on which the <TT>FEValues</TT> object was last reinitialized.
<P>
Given this, the innermost loop of <TT>assemble_system</TT> computes the
local contributions to the matrix in the following elegant way (the variable
<TT>stress_strain_tensor</TT>, corresponding to the tensor <IMG
- WIDTH="16" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img17.png"
- ALT="$ C$">, has
+ ALT="$ C$">
+, has
previously been initialized with the result of the first function above):
<PRE>
for (unsigned int i=0; i<dofs_per_cell; ++i)
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="331" HEIGHT="72" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="331" HEIGHT="71" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img88.png"
ALT="\begin{gather*}\begin{split}f^K_i &= (\vec f, \varphi_i)_K -(\sigma^{n-1},\varep...
...gma^{n-1}_q : \varepsilon(\varphi_i(\vec x_q)) \right\} \end{split}\end{gather*}"></TD>
<BR CLEAR="ALL"><P></P>
to the right hand side of (<A HREF="#eq:linear-system">17</A>) is equally
straightforward (note that we do not consider any boundary tractions <IMG
- WIDTH="13" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img89.png"
- ALT="$ \vec b$"> here). Remember that we only had to store the old stress in the
+ ALT="$ \vec b$">
+ here). Remember that we only had to store the old stress in the
quadrature points of cells. In the program, we will provide a variable
<TT>local_quadrature_points_data</TT> that allows to access the stress
<!-- MATH
$\sigma^{n-1}_q$
-->
<IMG
- WIDTH="38" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="39" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img90.png"
- ALT="$ \sigma^{n-1}_q$"> in each quadrature point. With this the code for the right
+ ALT="$ \sigma^{n-1}_q$">
+ in each quadrature point. With this the code for the right
hand side looks as this, again rather elegant:
<PRE>
for (unsigned int i=0; i<dofs_per_cell; ++i)
$\vec f(\vec x_q) \cdot \varphi_i(\vec x_q)$
-->
<IMG
- WIDTH="96" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img91.png"
- ALT="$ \vec f(\vec x_q) \cdot \varphi_i(\vec x_q)$">, we have made use of the fact that for the chosen finite element, only
+ ALT="$ \vec f(\vec x_q) \cdot \varphi_i(\vec x_q)$">
+, we have made use of the fact that for the chosen finite element, only
one vector component (namely <TT>component_i</TT>) of <IMG
- WIDTH="19" HEIGHT="29" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="20" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img83.png"
- ALT="$ \varphi_i$"> is
+ ALT="$ \varphi_i$">
+ is
nonzero, and that we therefore also have to consider only one component of
<!-- MATH
$\vec f(\vec x_q)$
-->
<IMG
- WIDTH="39" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="41" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img92.png"
- ALT="$ \vec f(\vec x_q)$">.
+ ALT="$ \vec f(\vec x_q)$">
+.
<P>
This essentially concludes the new material we present in this function. It
$\Delta \vec u^n$
-->
<IMG
- WIDTH="35" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
+ WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img39.png"
- ALT="$ \Delta \vec u^n$"> computed before, we update the stress values in all quadrature points
+ ALT="$ \Delta \vec u^n$">
+ computed before, we update the stress values in all quadrature points
according to (<A HREF="#eq:stress-update">18</A>) and (<A HREF="#eq:stress-update+rot">20</A>),
including the rotation of the coordinate system.
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="310" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="311" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img93.png"
ALT="$\displaystyle E(\vec u) = (\varepsilon(\vec u), C\varepsilon(\vec u))_{\Omega} - (\vec f, \vec u)_{\Omega} - (\vec b, \vec u)_{\Gamma_N},$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="86" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="88" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img94.png"
ALT="$\displaystyle f(\sigma(\vec u)) \le 0$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
- WIDTH="100" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="101" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img95.png"
ALT="$\displaystyle \nu = \frac{\lambda}{2(\lambda+\mu)},$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
$\lambda,\mu$
-->
<IMG
- WIDTH="29" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="30" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img96.png"
- ALT="$ \lambda,\mu$"> are the Lamé constants of the material.
+ ALT="$ \lambda,\mu$">
+ are the Lamé constants of the material.
Physical constraints indicate that <!-- MATH
$-1\le \nu\le \tfrac 12$
-->
<IMG
- WIDTH="85" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="86" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img97.png"
- ALT="$ -1\le \nu\le \tfrac 12$">. If <IMG
- WIDTH="12" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ ALT="$ -1\le \nu\le \tfrac 12$">
+. If <IMG
+ WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
SRC="step-18.data/intro/img98.png"
ALT="$ \nu$">
+
approaches <IMG
- WIDTH="13" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
+ WIDTH="15" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img99.png"
- ALT="$ \tfrac 12$">, then the material becomes incompressible. In that
+ ALT="$ \tfrac 12$">
+, then the material becomes incompressible. In that
case, pure displacement-based formulations are no longer appropriate for the
solution of such problems, and stabilization techniques have to be employed
for a stable and accurate solution. The book and paper cited above give
transfer the data that is stored in the quadrature points of the cells of the
old mesh to the new mesh, preferably by some sort of projection scheme. This
is only slightly messy in the sequential case; in fact, the functions
-<tt>FETools::get_projection_from_quadrature_points_matrix</tt> will do
-the projection, and the <tt>FiniteElement::get_restriction_matrix</tt> and
-<tt>FiniteElement::get_prolongation_matrix</tt> functions will do the
-transfer between mother and child cells. However, it becomes complicated once
-we run the program in parallel, since then each process only stores this data
-for the cells it owned on the old mesh, and it may need to know the values of
-the quadrature point data on other cells if the corresponding cells on the new
-mesh are assigned to this process after subdividing the new mesh. A global
-communication of these data elements is therefore necessary, making the entire
-process a little more unpleasant.
+<TT>FETools</TT> <TT>::</TT> <TT>get_projection_from_quadrature_points_matrix</TT> will do
+the projection, and the <TT>FiniteElement</TT> <TT>::</TT> <TT>get_restriction_matrix</TT> and
+<TT>FiniteElement</TT> <TT>::</TT> <TT>get_prolongation_matrix</TT> functions will do the
+transfer between mother and child cells. However, it becomes complicated
+once we run the program in parallel, since then each process only stores this
+data for the cells it owned on the old mesh, and it may need to know the
+values of the quadrature point data on other cells if the corresponding cells
+on the new mesh are assigned to this process after subdividing the new mesh. A
+global communication of these data elements is therefore necessary, making the
+entire process a little more unpleasant.
<P>
Conversely, if you want to run the program in 2d (after making the necessary
changes to accommodate for a 2d geometry), you have to change the Makefile
back to allow for 2d.
-
projects / dealii-svn.git / commitdiff
