OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
f0f4e399ab59f9788b840526f0af0ac7f47af0b7
IFEM/doc/Tutorials/GettingStarted.tex
Knut Morten Okstad f0f4e399ab Updated the ASM class hierarchy in the Getting Started presentation. 
Updated/corrected the Poisson tutorial.
2016-01-31 14:19:08 +01:00

1150 lines
37 KiB
TeX
Raw Blame History

% $Id$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass{beamer}
\usepackage[english]{babel}
\usepackage{eqalign}
\usepackage{bm}
\usepackage{epic}
\usepackage{epsfig}
\def\hideframe{\frame<article>}
\def\classitem#1{\item{\color{red}\tt#1} -}
\def\hugebox#1{\colorbox{white}{\fboxrule=0.5mm\fbox{\huge~#1~}}}
\def\BigBox#1#2{\fboxrule=0.3mm\fbox{\begin{minipage}{42mm}
{\small\color{red}#1} -- \scriptsize #2\end{minipage}}}
\title{\Large{\sl IFEM} - getting started}
\author{Knut Morten Okstad}
\institute{SINTEF ICT, Department of Applied Mathematics}
\hypersetup{colorlinks=true,urlcolor=blue}
\hypersetup{pdfauthor=Knut Morten Okstad}
\hypersetup{pdfsubject=Getting started with IFEM as simulation tool}
\hypersetup{pdfkeywords=NURBS Isogeometry FEA Object-orientation}
\begin{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frame %1
{
\frametitle{\footnotesize ICADA - Note 2011-001:}
\titlepage
}
\frame %2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{IFEM module overview}
\begin{picture}(300,150)
\thicklines
\put(140, 66){\includegraphics[width=1cm]{user}}
\put(149, 60){\tiny user}
\put(117, 95){\vector(-2, 1){30}}\put(107,100){\vector( 2,-1){30}}
\put( 0,140){\BigBox{GEOM}{Geometry Modeller
\begin{itemize}\tiny
\item Edit g2-files (GoTools) by hand
\item Using Rhino and small C++ tools
\item Python-based modeling tool
\end{itemize}}}
\uncover<2->{
\put(190, 95){\vector( 2, 1){30}}\put(200,100){\vector(-2,-1){30}}
\put(135,145){\vector( 1, 0){60}}
\put(200,140){\BigBox{GPM}{Grid/Property Modeller
\begin{itemize}\tiny
\item Automatic global node numbering for multi-patch models
\item Interactive assignment of property codes to topological entities
\end{itemize}}}}
\uncover<3->{
\put(190, 65){\vector( 2,-1){30}}\put(200, 60){\vector(-2, 1){30}}
\put(260,113){\vector( 0,-1){65}}
\put(200, 10){\BigBox{SIM}{Numerical Simulation
\begin{itemize}\tiny
\item Object-oriented framework for Isogeometric Finite Element Analysis
\item The user has to program his/her own application
\item A few sample applications are provided (Poisson, Linear elasticity)
\end{itemize}}}}
\uncover<4->{
\put(117, 65){\vector(-2,-1){30}}\put(107, 60){\vector( 2, 1){30}}
\put(195, 13){\vector(-1, 0){60}}
\put( 0, 10){\BigBox{RESV}{Result Visualization
\begin{itemize}\tiny
\item Currently, GLview Inova (Ceetron) \\ is used based on VTF file format
\item HDF5 file format is used for result storage and restart, from which
converters to any preferred visualization tool can be made
\end{itemize}}}}
\end{picture}
}
\frame %3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Major class hierarchies of the SIMulation environment}
\begin{enumerate}
\item<1-> {\color{red}NonLinSIM} - Nonlinear simulation driver
\begin{itemize}
\item Administers time/load step loop of the solution algorithm
\item Newton iteration loop,
\footnotesize convergence check, configuration update
\end{itemize}
\item<2-> {\color{red}SystemMatrix/Vector} - Linear algebra system level
\begin{itemize}
\item Interface to equation solvers (direct/iterative, serial/parallel)
\item Sub-classes for various linear equation solver packages
\end{itemize}
\item<3-> {\color{red}SIMbase} - System level drivers
\begin{itemize}
\item Administering an assembly of spline patches (blocks)
\item Sub-classes for problem-specific input and setup
\end{itemize}
\item<4-> {\color{red}ASMbase} - Block/patch level drivers
\begin{itemize}
\item Administers the element loop and numerical integration loop within
a block (spline patch)
\item Sub-classes depending on discretization
(Splines/NURBS, Lagrange, Spectral)
\item Uses {\sl GoTools} to evaluate basis functions at integration points
{\scriptsize\url{http://www.sintef.no/Projectweb/Geometry-Toolkits/GoTools}}
\end{itemize}
\item<5-> {\color{red}Integrand} - Integration point level
\begin{itemize}
\item Administers the problem-dependent calculations at an integration point
(interior and boundary integrals)
\item Problem-specific sub-classes
\end{itemize}
\end{enumerate}
\begin{picture}(0,0)
\color{red}
\uncover<6->{\put(55,97){\hugebox{Isogeometric level}}}
\uncover<7->{\put(55,25){\hugebox{User/Application level}}}
\end{picture}
}
\frame %4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{ASM class hierarchy - \em the Isogeometry level}
\begin{picture}(310,205)
\put(240,195){\vector(1,0){15}}
\put(258,193){\scriptsize 'is-a' relationship}
\dashline[50]{2}(240,185)(255,185)\put(250,185){\vector(1,0){5}}
\put(258,183){\scriptsize 'has-a' relationship}
\put( 90,190){\color{red}\fbox{\small ASMbase}}
\uncover<2->{
\put( 20,150){\color{red}\fbox{\small ASMstruct}}
\put( 47,161){\vector(2,1){50}}}
\uncover<2-8>{
\put(120,150){\color{green}\fbox{\small GeomObject}}
\dashline[50]{2}(72,154)(120,154)\put(115,154){\vector(1,0){5}}}
\uncover<3->{
\put( 37,121){\vector(1,3){8}}
\put( 12,110){\color{red}\fbox{\small ASMs1D}}}
\uncover<3-3>{
\put(110,110){\color{green}\fbox{\small SplineCurve}}
\put(140,120){\vector(1,3){8}}
\dashline[50]{2}(56,114)(110,114)\put(105,114){\vector(1,0){5}}}
\uncover<4->{\put( 15,105){\colorbox{white}{\color{red}\fbox{\small ASMs2D}}}}
\uncover<4-4>{
\put(110,108){\color{green}\fbox{\small SplineSurface}}
\put(140,119){\vector(1,3){8}}
\dashline[50]{2}(62,110)(110,112)\put(105,112){\vector(1,0){5}}}
\uncover<5->{\put(20,100){\colorbox{white}{\color{red}\fbox{\small ASMs3D}}}}
\uncover<5-8>{
\put(110,106){\color{green}\fbox{\small SplineVolume}}
\put(140,118){\vector(1,3){8}}
\dashline[50]{2}(67,105)(110,110)\put(105,110){\vector(1,0){5}}}
\uncover<6->{
\put( 43,71){\vector(0,1){25}}
\put( 15,60){\color{red}\fbox{\small ASMs3Dmx}}}
\uncover<6-8>{\dashline[50]{2}(72,65)(133,98)\put(130,95){\vector(1,1){6}}}
\uncover<7->{
\put(130,60){\vector(-2,1){72}}
\put(100,50){\color{red}\fbox{\small ASMs3DLag}}}
\uncover<8->{
\put( 70,20){\vector(2,1){48}}
\put( 30,10){\color{red}\fbox{\small ASMs3DmxLag}}}
\uncover<9->{
\put(210,150){\color{red}\fbox{\small ASMunstruct}}
\put(130,150){\color{green}\fbox{\small LRSpline}}
\put(235,161){\vector(-4,1){100}}
\dashline[50]{2}(210,153)(173,153)\put(178,153){\vector(-1,0){5}}}
\uncover<10->{
\put(250,121){\vector(-1,3){8}}
\put(150,121){\vector(0,1){23}}
\put(222,110){\color{red}\fbox{\small ASMu2D}}}
\uncover<10-10>{
\dashline[50]{2}(220,113)(200,113)\put(199,113){\vector(-1,0){5}}
\put(120,110){\color{green}\fbox{\small LRSplineSurface}}}
\uncover<11->{
\put(212,105){\colorbox{white}{\color{red}\fbox{\small ASMu3D}}}
\put(110,108){\color{green}\fbox{\small LRSplineVolume}}
\dashline[50]{2}(214,109)(188,111)\put(190,111){\vector(-1,0){5}}}
\uncover<12->{
\put(238,70){\vector(0,1){31}}
\put(210,60){\color{red}\fbox{\small ASMu3Dmx}}}
\put(260,10){\color{red}\scriptsize{\sl IFEM} classes}
\uncover<2-8>{
\put(260,0){\color{green}\scriptsize{\sl GoTools} classes}}
\uncover<9->{
\put(260,0){\color{green}\scriptsize{\sl LRSpline} classes}}
\end{picture}
\uncover<14->{\small
{\color{red} ASMbase} is the interface between the isogeometric FE procedures,
and the solution algorithms (from above) and the physical problem to be solved
(from below).}
}
\begin{frame}[fragile] %5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{The main ASM methods}
\tiny
\begin{verbatim}
class ASMbase
{
public:
//! \brief Evaluates an integral over the interior patch domain.
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time) = 0;
//! \brief Evaluates a boundary integral over a patch face/edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index of the boundary face/edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time) = 0;
};
\end{verbatim}
\footnotesize
{\color{red}LocalIntegral} and {\color{red}GlobalIntegral} are interfaces to
the element-level and system-level matrices of the FE problem.
{\color{red}TimeDomain} contains the integration parameters
needed for nonlinear and/or time-dependent simulations.
\end{frame}
\frame %5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Numerical integration method for a 3D spline patch}
\scriptsize\tt
\begin{tabbing} bo\=ol\=AS\=Ms\=3D::integrateXX(\= \kill
bool {\color{red}ASMs3D}::{\color{blue}integrate}
({\color{red}Integrand}\& integrand, {\color{red}GlobalIntegral}\& glInt, \\
\>\>\>\>\> const {\color{red}TimeDomain}\& time) \\
\{\\
\uncover<2->{\>
{\rm Compute parameter values (u,v,w)
of all integration points within the patch} \+\\[1mm]
basis->{\color{green}SplineVolume}::%
{\color{blue}computeBasisGrid}(u,v,w,splineData); \+\\[1pt]
splineData {\rm contains derivatives w.r.t. u,v,w of all basis functions} \\
{\rm at all integration points and the function values themselves}\-\-\\[1mm]}
\uncover<3->{
{\rm Loop over elements (knot-spans);} do iel=0,nel-1 \\}
\uncover<4->{\>\>
If {\rm current knot span is non-zero in all three directions} then \\}
\uncover<5->{\>\>\>
{\color{red}LocalIntegral}* A =
{\color{blue}integrand}.getLocalIntegral(iel); \\\>\>\>
{\rm Initialize for numerical integration over the element} \\\>\>\>
{\rm Fetch nodal coordinates (control points) for current element,} Xnod \\}
\uncover<6->{\>\>\>
{\rm Loop over integration points;}
do i=1,nGauss, j=1,nGauss, k=1,nGauss \\}
\uncover<7->{\>\>\>\>
{\rm Fetch data from} splineData
{\rm belonging to current integration point;} N, dN/du \\}
\uncover<8->{\>\>\>\>
{\rm Compute Cartesian coordinates and Jacobian;}
X = N*Xnod, J = dN/du*Xnod \\\>\>\>\>
{\rm and the gradient;} dN/dX = dN/du * J$^{-1}$ \\}
\uncover<9->{\>\>\>\>
integrand.{\color{blue}evalInt}%
(*A, time, detJ*weight, N, dN/dX, X); \\}
\uncover<6->{\>\>\>
end do i, j, k \\}
\uncover<10->{\>\>\>
integrand.{\color{blue}finalizeElement}(*A); \\}
\uncover<11->{\>\>\>
glInt.{\color{blue}assemble}(A->ref(), MGEL[iel]); \\\>\>\>
A->{\color{blue}destruct}(); \\}
\uncover<4->{\>\>
end if \\}
\uncover<3->{\>
end do iel \\}
\}
\end{tabbing}
\uncover<12->{
\thicklines\color{red}
\begin{picture}(0,0)
\put( 6,224){\line(1,0){290}}
\put( 6,182){\line(1,0){290}}
\put( 6,182){\line(0,1){42}}
\put(296,182){\line(0,1){42}}
\end{picture}}
}
\begin{frame}[fragile] %7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{The ``user'' interface ...}
\tiny
\begin{verbatim}
class Integrand
{
//! \brief Evaluates the integrand at an interior point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalInt(LocalIntegral& elmInt, const FiniteElement& fe,
const TimeDomain& time, const Vec3& X) const;
//! \brief Evaluates the integrand at a boundary point.
//! \param elmInt The local integral object to receive the contributions
//! \param[in] fe Finite element data of current integration point
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X Cartesian coordinates of current integration point
//! \param[in] normal Boundary normal vector at current integration point
//!
//! \details The default implementation forwards to the stationary version.
//! Reimplement this method for time-dependent or non-linear problems.
virtual bool evalBou(LocalIntegral& elmInt, const FiniteElement& fe,
const TimeDomain& time,
const Vec3& X, const Vec3& normal) const;
. . .
};
\end{verbatim}
\vskip-5mm\footnotesize
Overloaded versions of these method interfaces exist without the
{\tt TimeDomain} argument, for stationary/linear problems.
\end{frame}
\begin{frame}[fragile] %8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Finite element data at integration point level}
\scriptsize
\begin{verbatim}
class FiniteElement
{
public:
int iel; //!< Element identifier
size_t iGP; //!< Global integration point counter
double u; //!< First parameter of current point
double v; //!< Second parameter of current point
double w; //!< Third parameter of current point
double xi; //!< First local coordinate of current integration point
double eta; //!< Second local coordinate of current integration point
double zeta; //!< Third local coordinate of current integration point
double h; //!< Characteristic element size
Vector N; //!< Basis function values
Vector Navg; //!< Volume-averaged basis function values
Matrix dNdX; //!< First derivatives (gradient) of the basis functions
Matrix3D d2NdX2; //!< Second derivatives of the basis functions
Matrix G; //!< Matrix used for stabilized methods
double detJxW; //!< Weighted determinant of the coordinate mapping
};
\end{verbatim}
An object of this class is used to transport all integration point quantities
to the application-dependent integrands.
\end{frame}
\frame %6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Framework for two-field mixed formulations}
\begin{itemize}
\item<1-> Two sets of basis functions --
the first basis should be of one order higher than the second
\begin{itemize}
\item<2-> First approach: Establish the first basis by order-elevating
the second basis (yields only ${\cal C}^{p-2}$ continuity
for the first basis).
\item<2-> Second approach: Add one extra control point for the first basis
in each parameter direction, and then reparameterize (both bases
will posess ${\cal C}^{p-1}$ continuity but will have separate
control point locations).
\end{itemize}
\item<3-> The knot-span elements become the same for the two bases,
$\Rightarrow$ simplifies the finite element topology management.
\item<4-> Since the geometry represented by the two bases will be identical,
it suffice to use the second (lowest order) basis only,
when evaluating the Jacobian of the geometry mapping
and the basis function gradients w.r.t.\ Cartesian coordinates.
\item<5-> The user only needs to relate to the lowest-order grid/basis, the
higher order basis is established internally automatically.
\end{itemize}
}
\frame %7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Numerical integration for the two-field mixed method}
\scriptsize\tt
\begin{tabbing} bo\=ol\=AS\=Ms\=3Dmx::integrateXX(\= \kill
bool {\color{red}ASMs3Dmx}::{\color{blue}integrate}
({\color{red}Integrand}\& integrand, {\color{red}GlobalIntegral}\& glInt, \\
\>\>\>\>\> const {\color{red}TimeDomain}\& time) \\
\{\\
\uncover<2->{\>
{\rm Compute parameter values (u,v,w)
of all integration points within the patch} \+\\[1mm]
basis1->{\color{green}SplineVolume}::%
{\color{blue}computeBasisGrid}(u,v,w,splineData1); \\
basis2->{\color{green}SplineVolume}::%
{\color{blue}computeBasisGrid}(u,v,w,splineData2); \-\\[1mm]}
\uncover<3->{\>
{\rm Loop over elements (knot-spans);} do iel=0,nel-1 \+\+\\
If {\rm current knot span is non-zero in all three directions} then \-\-\\}
\uncover<4->{\>\>\>
{\color{red}LocalIntegral}* A =
{\color{blue}integrand}.getLocalIntegral(iel); \\\>\>\>
{\rm Initialize for numerical integration over the element} \\\>\>\>
{\rm Fetch nodal coordinates for current element,} Xnod
{\rm (Note: for }basis2 {\rm only)} \\}
\uncover<5->{\>\>\>
{\rm Loop over integration points;}
do i=1,nGauss, j=1,nGauss, k=1,nGauss \\}
\uncover<6->{\>\>\>\>
{\rm Fetch data from} splineData[12]; N$_1$, dN$_1$/du, N$_2$, dN$_2$/du\\}
\uncover<7->{\>\>\>\>
{\rm Compute Cartesian coordinates and Jacobian;}
X = N$_2$*Xnod, J = dN$_2$/du*Xnod \\\>\>\>\>
{\rm and the gradients;} dN$_1$/dX = dN$_1$/du * J$^{-1}$,
dN$_2$/dX = dN$_2$/du * J$^{-1}$, \\}
\uncover<8->{\>\>\>\>
integrand.{\color{blue}evalInt}(*A, time, detJ*weight, \\\>\>\>\>\>
N$_1$, dN$_1$/dX, N$_2$, dN$_2$/dX, X); \\}
\uncover<5->{\>\>\>
end do i, j, k \\}
\uncover<9->{\>\>\>
integrand.{\color{blue}finalizeElement}(*A); \\\>\>\>
glInt.{\color{blue}assemble}(A->ref(), MGEL[iel]); \\\>\>\>
A->{\color{blue}destruct}(); \\}
\uncover<3->{\>\>
end if \\\>
end do iel \\}
\}
\end{tabbing}
}
\frame %8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Current Applications - {\color{red}\tt Integrand} sub-classes}
\small
\begin{itemize}
\classitem{Poisson} Simple scalar equation
\pause
\classitem{Elasticity} Solid mechanics problems
\begin{itemize}
\classitem{LinearElasticity} Linear elasticity, isotropic material
\classitem{NonlinearElasticityTL} Finite deformation elasticity, \\
Total Lagrangian formulation, linear elastic material
\classitem{NonlinearElasticityUL} Finite deformation elasticity,
Updated Lagrangian formulation, linear elastic, Neo-Hookean and
plastic materials (linear elastic only in public version)
\pause
\begin{itemize}
\classitem{NonlinearElasticityULMX} Incompressible and nearly incompressible
materials, mixed formulation with internal pressure and volumetric change
modes (not in public version)
\classitem{NonlinearElasticityULMixed} ...,
mixed formulation with continuous pressure and volumetric change fields
(not public) \\
\classitem{NonlinearElasticityFbar} ...,
$\bar{F}$-formulation for the handling of nearly-incompressible materials
(not public)
\end{itemize}
\end{itemize}
\pause
\item Navier--Stokes CFD solvers (not part of the {\sl ICADA\/} project)
\pause
\item Many others (coupled simulators, etc.)
\end{itemize}
}
\frame %9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Implementational issues (integration point level)}
\begin{itemize}
\item<1-> Using splines as basis function, especially the higher-order ones,
the ``elements'' become large (in terms of nodal connectivities)
$\Rightarrow$ large, dense element matrices
\item<2-> Element-level linear algebra: Use machine-optimized {\bf BLAS}
rather than inline C++ code
\item<2-> Important to express the nonlinear FE formulation on {\em matrix}
form (Voigt notation) --- not {\em tensor} form
\item<3-> In addition: Integration and assembly of element-level matrices
is done in parallel using multi-threading (OpenMP)
\end{itemize}
}
\frame %10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{System-level linear algebra -- equation solving}
\begin{itemize}
\item Interfaced through classes {\color{red}\tt SystemMatrix} and
{\color{red}\tt SystemVector} with sub-classes for particular solvers.
\vskip5pt
\item Current available linear equation solvers:
\begin{itemize}
\item LAPACK {\tt DGESV} (dense matrices, small problems only)
\item SuperLU (direct methods)
\hfill{\scriptsize\url{http://crd.lbl.gov/~xiaoye/SuperLU}}
\item PETSc (iterative methods)
\hfill{\scriptsize\url{http://www.mcs.anl.gov/petsc}}
\item Parallelization on distributed memory based on PETSc/MPI
\end{itemize}
\end{itemize}
}
\frame %11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\frametitle{Detailed source code documentation}
See the doxygen-generated html-pages \url{../html/index.html}
}
\frame %12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{
\begin{center}
\huge\color{blue} Tutorial: Poisson equation in $R^2$
\end{center}
Given a heat source function $f(x,y)$ defined over a domain
$\Omega\in R^2$, a flux function $h(x,y)$ defined over the boundary
$\partial\Omega_h$, and a function $g(x,y)$ defined over the boundary
$\partial\Omega_g = \partial\Omega\setminus\partial\Omega_h$,
find the scalar function $u(x,y)$ satisfying
\begin{eqnarray}
\left.\eqalign{
q_{i,i} &= f \cr
q_i & = -\kappa_{ij} u_{,j}
}\right\} &\forall& \{x,y\}\in\overline{\Omega} \\[2mm]
q_i n_i = h &\forall& \{x,y\}\in\partial\Omega_h \\
u = g &\forall&\{x,y\}\in\partial\Omega_g
\end{eqnarray}
where $\kappa_{ij}$ is the conductivity tensor and $n_i$ defines the outward-
directed unit normal vector on $\partial\Omega_h$.
}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.58\textwidth}
\tiny\begin{verbatim}
class Poisson : public IntegrandBase
{
protected:
// Physical properties
double kappa; //!< Conductivity (constant)
RealFunc* fluxFld; //!< Boundary normal flux field
RealFunc* heatSrc; //!< Interior heat source field
\end{verbatim}
\column{0.42\textwidth} \small
Define the class {\tt Poisson} as an {\tt IntegrandBase} subclass
{\scriptsize(the class {\tt IntegrandBase} inherits {\tt Integrand})},
containing data and methods specific to the 2D Poisson problem
(assuming constant conductivity).
\end{columns}
\pause
\begin{columns}[c]
\column{0.58\textwidth}
\tiny\begin{verbatim}
public:
Poisson() : kappa(1.0), fluxFld(0), heatSrc(0)
{
primsol.resize(1);
}
virtual ~Poisson() {}
\end{verbatim}
\column{0.42\textwidth}\small
Class constructor and destructor.
The constructor {\tt Poisson()} initializes the data members.
\end{columns}
\pause
\begin{columns}[c]
\column{0.58\textwidth}
\tiny\begin{verbatim}
void setMaterial(double K) { kappa = K; }
void setFlux(RealFunc* ff) { fluxFld = ff; }
void setSource(RealFunc* src) { heatSrc = src; }
\end{verbatim}
\column{0.42\textwidth} \small
Initialization of physical properties.
\end{columns}
\pause
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
virtual LocalIntegral* getLocalIntegral(size_t nen, size_t,
bool neumann) const;
\end{verbatim}
\column{0.4\textwidth} \small
Virtual method returning an element matrix object for the Poisson integrand.
\end{columns}
\end{frame}
\begin{frame}[fragile]
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
virtual bool evalInt(LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X) const;
virtual bool evalBou(LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X,
const Vec3& normal) const;
virtual bool evalSol(Vector& s,
const FiniteElement& fe,
const Vec3& X,
const std::vector<int>& MNPC) const;
virtual bool evalSol(Vector& s,
const VecFunc& asol,
const Vec3& X) const;
\end{verbatim}
\column{0.4\textwidth} \small
Virtual methods for integrand and solution field evaluation.
\end{columns}
\pause
\vskip-5mm
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
virtual NormBase* getNormIntegrand(AnaSol* asol = 0) const;
bool evalSol(Vector& s,
const Vector& eV,
const Matrix& dNdX,
const Vec3& X) const;
bool formCmatrix(Matrix& C, const Vec3& X,
bool invers = false) const;
};
\end{verbatim}
\column{0.4\textwidth}\small
Methods for solution norm integration.
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
LocalIntegral* Poisson::getLocalIntegral (size_t nen, size_t,
bool neumann) const
{
ElmMats* result = new ElmMats();
result->rhsOnly = neumann;
result->withLHS = !neumann;
result->resize(neumann ? 0 : 1, 1);
result->redim(nen);
return result;
}
\end{verbatim}
\column{0.4\textwidth}\small
Element initialization:
Allocate an element matrix object and set the size of the matrices
based on the number of element nodes.
Indicate whether the left-hand-side matrices are to be integrated or not.
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
bool Poisson::evalInt (LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X) const
{
ElmMats& elMat = static_cast<ElmMats&>(elmInt);
if (!elMat.A.empty())
{
Matrix C; C.diag(kappa,2); // Diagonal constitutive matrix
Matrix CB;
CB.multiply(C,fe.dNdX,false,true).multiply(fe.detJxW);
elMat.A.front().multiply(fe.dNdX,CB,false,false,true);
}
if (heatSrc && !elMat.b.empty())
elMat.b.front().add(fe.N,(*heatSrc)(X)*fe.detJxW);
return true;
}
bool Poisson::evalBou (LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X, const Vec3& normal) const
{
ElmMats& elMat = static_cast<ElmMats&>(elmInt);
if (!fluxFld || elMat.b.empty()) return false;
double h = -(*fluxFld)(X); // normal flux at point X
elMat.b.front().add(fe.N,h*fe.detJxW);
return true;
}
\end{verbatim}
\column{0.4\textwidth}\small
Integrand evaluations:
Assuming here $\kappa_{ij} = \kappa\delta_{ij}$
$$
[CB] = [C]\cdot[\partial N/\partial\bm{X}]^T|J|w
$$
$$
[eM] = \sum\left([\partial N/\partial\bm{X}]\cdot[CB]\right)
$$
$$
\{eS\} = \sum\left(f\{N\}|J|w\right)
$$
\vskip15mm
$$
\{eS\} = \sum\left(-h\{N\}|J|w \right)
$$
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.55\textwidth}
\tiny\begin{verbatim}
bool Poisson::evalSol (Vector& q,
const FiniteElement& fe, const Vec3& X,
const std::vector<int>& MNPC) const
{
if (primsol.front().empty()) return false;
Vector eV;
int ierr = utl::gather(MNPC,1,primsol.front(),eV);
if (ierr > 0) return false;
Matrix C; C.diag(kappa,2); // Diagonal constitutive matrix
// Evaluate the heat flux vector
Matrix CB;
CB.multiply(C,fe.dNdX,false,true).multiply(eV,q);
q *= -1.0;
return true;
}
bool Poisson::evalSol (Vector& s, const VecFunc& asol,
const Vec3& X) const
{
s = Vector(asol(X).ptr(),2);
return true;
}
\end{verbatim}
\column{0.45\textwidth}\small
Secondary solution evaluation:
$$
\bm{q} = [C]\cdot[\partial N/\partial\bm{X}]^T\cdot\{eV\}
$$
\vskip3cm
Analytical solution
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
class PoissonNorm : public NormBase
{
VecFunc* anasol; //!< Analytical heat flux
public:
PoissonNorm(Poisson& p, VecFunc* a = 0)
: NormBase(p), anasol(a) {}
virtual ~PoissonNorm() {}
virtual bool hasBoundaryTerms() const { return true; }
virtual bool evalInt(LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X) const;
virtual bool evalBou(LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X,
const Vec3& normal) const;
};
NormBase* Poisson::getNormIntegrand (AnaSol* asol) const
{
if (asol)
return new PoissonNorm(*const_cast<Poisson*>(this),
asol->getScalarSecSol());
else
return new PoissonNorm(*const_cast<Poisson*>(this));
}
\end{verbatim}
\column{0.45\textwidth}\small
Accompanying class for solution norm integration \vskip\baselineskip
{\tt NormBase} is a sub-class of {\tt Integrand} with a couple of added
methods common to all norm classes.
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
bool PoissonNorm::evalInt (LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X) const
{
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = static_cast<ElmNorm&>(elmInt);
// Evaluate the inverse constitutive matrix at this point
Matrix Cinv;
problem.formCmatrix(Cinv,X,true);
// Evaluate the finite element heat flux field
Vector sigmah;
problem.evalSol(sigmah,pnorm.vec.front(),fe.dNdX,X);
// Integrate the energy norm a(u^h,u^h)
pnorm[0] += sigmah.dot(Cinv*sigmah)*fe.detJxW;
// Integrate the external energy (h,u^h)
double u = pnorm.vec.front().dot(fe.N);
pnorm[1] += problem.getHeat(X)*u*fe.detJxW;
if (anasol) {
// Evaluate the analytical heat flux
Vector sigma((*anasol)(X).ptr(),sigmah.size());
// Integrate the energy norm a(u,u)
pnorm[2] += sigma.dot(Cinv*sigma)*fe.detJxW;
// Integrate the error in energy norm a(u-u^h,u-u^h)
sigma -= sigmah;
pnorm[3] += sigma.dot(Cinv*sigma)*fe.detJxW;
}
return true;
}
\end{verbatim}
\column{0.45\textwidth}\small
Norm integrand evaluation
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
bool PoissonNorm::evalBou (LocalIntegral& elmInt,
const FiniteElement& fe,
const Vec3& X,
const Vec3& normal) const
{
Poisson& problem = static_cast<Poisson&>(myProblem);
ElmNorm& pnorm = static_cast<ElmNorm&>(elmInt);
// Evaluate the surface heat flux
double t = problem.getTraction(X,normal);
// Evaluate the temperature field
double u = pnorm.vec.front().dot(fe.N);
// Integrate the external energy (t,u^h)
pnorm[1] += t*u*fe.detJxW;
return true;
}
\end{verbatim}
\column{0.45\textwidth}\small
Integration of external energy due to boundary flux.
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.55\textwidth}
\tiny\begin{verbatim}
class SIMPoisson2D : public SIM2D
{
Poisson prob; //!< Poisson data and methods
RealArray mVec; //!< Material data
public:
SIMPoisson2D() : SIM2D(1) { myProblem = &prob; }
virtual ~SIMPoisson2D() { myProblem = 0; }
protected:
virtual bool parse(const TiXmlElement* elem);
virtual bool initMaterial(size_t propInd);
virtual bool initNeumann(size_t propInd);
};
bool SIMPoisson2D::initMaterial (size_t propInd)
{
if (propInd >= mVec.size()) return false;
prob.setMaterial(mVec[propInd]);
return true;
}
bool SIMPoisson2D::initNeumann (size_t propInd)
{
SclFuncMap::const_iterator sit = myScalars.find(propInd);
if (sit == myVectors.end()) return false;
prob.setFlux(tit->second);
return true;
}
\end{verbatim}
\column{0.5\textwidth}\small
Simulation driver for 2D problems.\\[5mm]
\pause
Alternative, use templates to support multiple dimensions:
\tiny\begin{verbatim}
template<class Dim>
class SIMPoisson : public Dim
{
...
};
\end{verbatim}
\small where {\tt Dim} can be either {\tt SIM1D}, {\tt SIM2D} or {\tt SIM3D}.
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\tiny\begin{verbatim}
bool SIMPoisson2D::parse (const TiXmlElement* elem)
{
if (strcasecmp(elem->Value(),"poisson"))
return this->SIM2D::parse(elem);
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if (!strcasecmp(child->Value(),"isotropic")) {
int code = this->parseMaterialSet(child,mVec.size());
double kappa = 1000.0;
utl::getAttribute(child,"kappa",kappa);
if (code == 0)
prob.setMaterial(kappa);
mVec.push_back(kappa);
}
\end{verbatim}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\tiny\begin{verbatim}
else if (!strcasecmp(child->Value(),"source")) {
int code = -1; // Reserve negative code(s) for the source term function
while (myScalars.find(code) != myScalars.end()) --code;
std::string type;
utl::getAttribute(child,"type",type,true);
if (type == "expression" && child->FirstChild()) {
std::cout <<"\tHeat source function: "
<< child->FirstChild()->Value() << std::endl;
myScalars[code] = new EvalFunction(child->FirstChild()->Value());
}
else {
std::cerr <<" ** SIMPoisson2D::parse: Invalid source function "<< type << std::endl;
continue;
}
prob.setSource(myScalars[code]);
}
else if (!strcasecmp(child->Value(),"anasol")) {
std::string type;
utl::getAttribute(child,"type",type,true);
if (type == "expression") {
std::cout <<"\tAnalytical solution: Expression"<< std::endl;
if (!mySol) mySol = new AnaSol(child);
}
else
std::cerr <<" ** SIMPoisson2D::parse: Invalid analytical solution "<< type << std::endl;
}
return true;
}
\end{verbatim}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation}
\begin{columns}[c]
\column{0.6\textwidth}
\tiny\begin{verbatim}
int main (int argc, char** argv)
{
// (Lots of initialisations skipped here...)
// Read in model definitions and establish the FE data structures
SIMbase* model = new SIMPoisson2D(); // (or new SIMPoisson<SIM2D>;)
if (!model->read(infile))
return 1;
if (!model->preprocess(ignoredPatches,fixDup))
return 1;
model->setQuadratureRule(nGauss);
Matrix eNorm;
Vector gNorm, sol;
model->initSystem(solver,1,1);
model->setAssociatedRHS(0,0);
if (!model->assembleSystem())
return 2;
// Solve the linear system of equations
if (!model->solveSystem(sol,1))
return 3;
// Evaluate solution norms
if (!model->solutionNorms(Vectors(1,sol),eNorm,gNorm))
return 4;
// Print output to terminal and VTF, etc.
}
\end{verbatim}
\column{0.45\textwidth}\small
Core parts of the main program
\end{columns}
\end{frame}
\begin{frame}[fragile] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\frametitle{Tutorial: Poisson equation, sample input files}
\begin{columns}[c]
\column{0.7\textwidth}
\tiny\begin{verbatim}
<simulation>
<!-- General - geometry definitions !-->
<geometry>
<patchfile>square2D.g2</patchfile>
<raiseorder patch="1" u="2" v="2"/>
<refine type="uniform" patch="1" u="7" v="7"/>
<topologysets>
<set name="Dirichlet" type="edge">
<item patch="1">4</item>
</set>
<set name="Neumann" type="edge">
<item patch="1">3</item>
</set>
</topologysets>
</geometry>
<!-- General - boundary conditions !-->
<boundaryconditions>
<dirichlet set="Dirichlet" comp="1"/>
<neumann type="anasol" set="Neumann" comp="1"/>
</boundaryconditions>
<!-- Problem-specific block !-->
<poisson>
<source type="expression">PI*PI*cos(PI*x)*(2-y)</source>
<anasol type="expression">
<primary>cos(PI*x)*(2-y)</primary>
<secondary>PI*sin(PI*x)*(2-y)|cos(PI*x)</secondary>
</anasol>
</poisson>
</simulation>
\end{verbatim}
\column{0.3\textwidth}
\tiny square2d.g2:
\begin{verbatim}
200 1 0 0
3 0
2 2
0 0 1 1
2 2
0 0 1 1
0.0 0.0 0.0
2.0 0.0 0.0
0.0 2.0 0.0
2.0 2.0 0.0
\end{verbatim}
\pause
\small This is equivalent to both of the following:
\tiny\begin{verbatim}
<geometry scale="2.0"/>
<geometry Lx="2.0" Ly="2.0"/>
\end{verbatim}
\small and then {\scriptsize\tt <patchfile>} is not needed.
\end{columns}
\end{frame}
\end{document} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Reference in New Issue View Git Blame Copy Permalink