# you switch from debug to optimized mode). To avoid messing with
# -Wl,-rpath in our linker flags, we create one convenience library
# for all PETSc libraries in our own lib directory, and have this one
-# twiddle with whatever is necessary to link in PETSc.
+# twiddle with whatever is necessary to link in PETSc. This is only
+# used for PETSc <=3.0. After this version, PETSc creates one library
+# we can link to directly. This avoids problems with recent versions
+# of gcc.
ifeq ($(USE_CONTRIB_PETSC),yes)
ifeq ($(DEAL_II_PETSC_VERSION_MAJOR),2)
lib-contrib-petsc-path.g = $(DEAL_II_PETSC_DIR)/lib/$(DEAL_II_PETSC_ARCH)
lib-contrib-petsc-path.o = $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib
endif
- # First up is a twiddle for Darwin/MacOSX
- ifeq ($(findstring darwin,$(TARGET)),)
+ ifeq ($(enable-shared),yes)
+ # Starting with PETSc versions 2.3.*:
+ ifeq ($(DEAL_II_PETSC_VERSION_MAJOR),2)
+ lib-contrib-petsc.g = $(shell echo $(DEAL_II_PETSC_DIR)/lib/$(DEAL_II_PETSC_ARCH)/*$(lib-suffix))
+ lib-contrib-petsc.o = $(shell echo $(DEAL_II_PETSC_DIR)/lib/$(DEAL_II_PETSC_ARCH)/*$(lib-suffix))
+ else
+ # which is the similar as for PETSc 3.0.0:
+ ifeq ($(DEAL_II_PETSC_VERSION_MAJOR)$(DEAL_II_PETSC_VERSION_MINOR),30)
+ lib-contrib-petsc.g = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/*$(lib-suffix))
+ lib-contrib-petsc.o = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/*$(lib-suffix))
+ else
+ # but after that (petsc-3.1++), we can use the simpler PETSc
+ # default "--with-single-library=1" like this:
+ lib-contrib-petsc.g = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/libpetsc$(lib-suffix))
+ lib-contrib-petsc.o = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/libpetsc$(lib-suffix))
+ endif
+ endif # if PETSC_VERSION
+ else
+ # and finally this goes for static libraries
lib-contrib-petsc.g = $(LIBDIR)/libpetscall.g$(lib-suffix)
lib-contrib-petsc.o = $(LIBDIR)/libpetscall$(lib-suffix)
- else
- ifeq ($(enable-shared),yes)
- # Starting with PETSc versions 2.3.*:
- ifeq ($(DEAL_II_PETSC_VERSION_MAJOR),2)
- lib-contrib-petsc.g = $(shell echo $(DEAL_II_PETSC_DIR)/lib/$(DEAL_II_PETSC_ARCH)/*$(lib-suffix))
- lib-contrib-petsc.o = $(shell echo $(DEAL_II_PETSC_DIR)/lib/$(DEAL_II_PETSC_ARCH)/*$(lib-suffix))
- else
- # which is the similar as for PETSc 3.0.0:
- ifeq ($(DEAL_II_PETSC_VERSION_MAJOR)$(DEAL_II_PETSC_VERSION_MINOR),30)
- lib-contrib-petsc.g = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/*$(lib-suffix))
- lib-contrib-petsc.o = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/*$(lib-suffix))
- else
- # but after that (petsc-3.1++), we can use the simpler PETSc
- # default "--with-single-library=1" like this:
- lib-contrib-petsc.g = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/libpetsc$(lib-suffix))
- lib-contrib-petsc.o = $(shell echo $(DEAL_II_PETSC_DIR)/$(DEAL_II_PETSC_ARCH)/lib/libpetsc$(lib-suffix))
- endif
- endif # if PETSC_VERSION
- else
- # and finally this goes for static libraries
- lib-contrib-petsc.g = $(LIBDIR)/libpetscall.g$(lib-suffix)
- lib-contrib-petsc.o = $(LIBDIR)/libpetscall$(lib-suffix)
- endif
endif
endif # USE_CONTRIB_PETSC
</p>
<ol>
+<li> Changed/fixed: Several operations on Trilinos vectors that change the
+elements of these vectors were allowed accidentally for vectors that have
+ghost elements. This is a source of errors because a change on one
+MPI process will not show up on a different processor. Consequently, we
+intended to disallow all functions that modify vectors with ghost elements
+but this was not enforced for all of these functions. This is now fixed
+but it may lead to errors if your code relied on the existing behavior. The
+way to work around this is to only ever modify fully distributed vectors,
+and then copy it into a vector with ghost elements.
+<br>
+(Wolfgang Bangerth, 2012/08/06)
+
<li> Changed: The argument material_id of the estimate-functions of
the KellyErrorEstimator class is now of type types::material_id_type
with default value types::invalid_material_id, instead of type
<ol>
+<li>
+New: officially added support for clang 3.1 or newer.
+<br>
+(Timo Heister and Wolfgang Bangerth, 2012/08/07)
+
+<li>
+Changed: PETSc linking now prefers to use the libpetsc.so generated
+by PETSc starting from version 3.1+. This fixes the problem
+of linker errors on recent gcc/ubuntu versions.
+<br>
+(Timo Heister, 2012/08/07)
+
<li>
Fixed: On some 64-bit systems, we build deal.II with the <code>-m64</code>
flag but forgot to build UMFPACK with this flag as well, leading to
Preprint, University of Heidelberg, Germany, 2011.
</li>
+ <li>
+ <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann/"
+ target="_top">R. Hartmann</a>
+ <br>
+ <strong>Higher-order and adaptive discontinuous Galerkin
+ methods with shock-capturing applied to transonic turbulent
+ delta wing flow
+ </strong>
+ <br>
+ 50<sup>th</sup> AIAA Aerospace Sciences Meeting, AIAA 2012-459, 2012.
+ </li>
+
+ <li>
+ <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann/"
+ target="_top">R. Hartmann</a>
+ <br>
+ <strong>Higher-order and adaptive Discontinuous Galerkin
+ methods applied to turbulent delta wing flow </strong>
+ <br>
+ In H. P. Kreplin, editor, New Results in Numerical and
+ Experimental Fluid Mechanics VIII: Contributions to the 17th
+ STAB/DGLR Symposium Berlin, Germany, 2010, Notes on
+ Numerical Fluid Mechanics and Multidisciplinary Design,
+ Springer. In press.
+ </li>
+
<li> T. Heister, G. Rapin
<br>
<strong>Efficient augmented Lagrangian-type preconditioning
pp. 3186-3197, 2012.
</li>
+ <li> M. Wallraff, T. Leicht, M. Lange-Hegermann
+ <br>
+ <strong>Numerical flux functions for Reynolds-averaged Navier-Stokes
+ and kω turbulence model computations with a line-preconditioned
+ p-multigrid discontinuous Galerkin solver</strong>
+ <br>
+ Int. J. Numer. Meth. Fluids, published online, 2012.
+ </li>
+
+ <li> M. Wallraff, T. Leicht, M. Lange-Hegermann
+ <br>
+ <strong>Numerical flux functions for RANS-kω computations with a
+ line-implicit p-multigrid DG solver</strong>
+ <br>
+ 50<sup>th</sup> AIAA Aerospace Sciences Meeting, AIAA 2012-458, 2012.
+ </li>
+
<li> J. Willems
<br>
<strong>Robust Multilevel Methods for General Symmetric Positive Definite Operators
SIAM J. of Numerical Analysis, vol. 49, pp. 917-944, 2011.
</li>
- <li> R. Hartmann
- <br>
- <strong>Higher-order and adaptive Discontinuous Galerkin
- methods applied to turbulent delta wing flow </strong>
- <br>
- In H. P. Kreplin, editor, New Results in Numerical and
- Experimental Fluid Mechanics VIII: Contributions to the 17th
- STAB/DGLR Symposium Berlin, Germany, 2010, Notes on
- Numerical Fluid Mechanics and Multidisciplinary Design,
- Springer. To appear.
- </li>
-
- <li>
- <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann/"
- target="_top">R. Hartmann</a>
- <br>
- <strong>Higher-order and adaptive discontinuous Galerkin
- methods with shock-capturing applied to transonic turbulent
- delta wing flow
- </strong>
- <br>
- 50thAIAA Aerospace Sciences Meeting, 2012. Submitted.
- </li>
-
<li> <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann"
target="_top">R. Hartmann</a>, J. Held, T. Leicht
<br>
J. Comput. Phys., vol. 230, pp. 4268-4284, 2011.
</li>
- <li> <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann"
- target="_top">R. Hartmann</a>, J. Held, T. Leicht, F. Prill
- <br>
- <strong>Error estimation and adaptive mesh refinement for aerodynamic flows
- </strong>
- <br>
- In N. Kroll, H. Bieler, H. Deconinck, V. Couallier, H. van der Ven and K. Sorensen, editors,
- ADIGMA - A European Initiative on the Development of Adaptive
- Higher-Order Variational Methods for Aerospace Applications,
- Springer, to appear, 2010.
- </li>
-
<li> T. Heister
<br>
<strong>A Massively Parallel Finite Element Framework with
refinement for aerodynamic flows.
</strong>
<br>
- 49<sup>th</sup>49th AIAA Aerospace Sciences Meeting, 2011.
+ 49<sup>th</sup> AIAA Aerospace Sciences Meeting AIAA 2011-212, 2011.
</li>
<li>
Aerospace Science and Technology, vol. 14, pp. 512-529, 2010.
</li>
+ <li> <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann"
+ target="_top">R. Hartmann</a>, J. Held, T. Leicht, F. Prill
+ <br>
+ <strong>Error estimation and adaptive mesh refinement for aerodynamic flows
+ </strong>
+ <br>
+ In N. Kroll, H. Bieler, H. Deconinck, V. Couallier, H. van der Ven
+ and K. Sorensen, editors, ADIGMA - A European Initiative on the
+ Development of Adaptive Higher-Order Variational Methods for
+ Aerospace Applications,
+ <br>
+ vol. 113 of Notes on numerical fluid mechanics and
+ multidisciplinary design, pp. 339-353, Springer, 2010.
+ </li>
+
<li> <a href="http://ganymed.iwr.uni-heidelberg.de/~hartmann"
target="_top">R. Hartmann</a>, P. Houston
<br>
developed and tested on the following platforms:
<ul>
<li>Linux on x86, with gcc 3.3, 3.4, 4.0, 4.1, 4.2, 4.3, 4.4 and 4.5; and
- Intel ICC 10.0, 11.0 and 12;
+ Intel ICC 10.0, 11.0 and 12; clang 3.1
<li>Linux on x86_64 (Opteron, AMD64, and Intel EM64T chips) using gcc
- 3.3, 3.4, 4.0, 4.1, 4.3, 4.4 and 4.5;
+ 3.3, 3.4, 4.0, 4.1, 4.3, 4.4 and 4.5; clang 3.1
<li>Mac OS X 10.2 or later, with gcc 4.0 or later (Please see the <a href="http://dealii.sourceforge.net/index.php?title=Main_Page"
target="_top">deal.II Wiki</a> for installation issues);
<li>Sun Solaris 2.5, 2.6, 7, 8, and 9 on Sparc, with
solution vector locally. The goal is to run this code on hundreds or maybe
even thousands of processors, at reasonable scalability.
+@note Even though it has a larger number, step-40 comes logically before the
+current program. You will probably want to look at step-40 before you try to
+understand the what we do here.
+
MPI is a rather awkward interface to program with. It is a semi-object
oriented set of functions, and while one uses it to send data around a
network, one needs to explicitly describe the data types because the MPI
--- /dev/null
+\documentclass{article}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+
+\title{Documentation of step-42, An obstacle problem for elasto-plastic material behavior in three dimensions}
+\author{Joerg Frohne}
+\date{Juni, 2012}
+
+\begin{document}
+
+\section{Introduction}
+
+This example composes an advanced version of step-41 since it considers an
+elasto-plastic material behavior with isotropic hardening in three dimensions.
+That means that we have to take care of an additional nonlinearity by the
+material behavior. An other difference compared to step-41 is that
+the contact area is arranged at the boundary of the deformable body now.\\
+Since you can slightly reach a few million degrees of freedom in three dimensions,
+even with adaptive mesh refinement, we decided to use trilinos and p4est to run
+our code in parallel. On the other hand we have to deal with hanging nodes because
+of the adaptive mesh, which is an other advance of step-41.
+
+
+\section{Classical formulation}
+
+The classical formulation of the problem possesses the following form:
+\begin{align*}
+ \varepsilon(u) &= A\sigma + \lambda & &\quad\text{in } \Omega,\\
+ \lambda(\tau - \sigma) &\geq 0\quad\forall\tau\text{ mit }\mathcal{F}(\tau)\leq 0 & &\quad\text{in } \Omega,\\
+ -\textrm{div}\ \sigma &= f & &\quad\text{in } \Omega,\\
+ u(\mathbf x) &= 0 & &\quad\text{on }\Gamma_D,\\
+ \sigma_t(u) &= 0,\quad\sigma_n(u)\leq 0 & &\quad\text{on }\Gamma_C,\\
+\sigma_n(u)(u_n - g) &= 0,\quad u_n(\mathbf x) - g(\mathbf x) \leq 0 & &\quad\text{on } \Gamma_C
+\end{align*}
+with $u\in H^2(\Omega)$. The vector valued function $u$ denotes the
+displacement in the deformable body. The first two lines describe the elast-plastic
+material behaviour. Therein the equation shows the deformation $\varepsilon (u)$ as the additive
+decomposition of the elastic part $A\sigma$ and the plastic part $\lambda$. $A$ is defined as
+the compliance tensor of fourth order which contains some material constants and $\sigma$ as the
+symmetric stress tensor of second order. So we have to consider the inequality in the second
+row component-by-component and furthermore we have to distinguish two cases.\\
+The continuous and convex function $\mathcal{F}$ denotes the von mises flow function
+$$\mathcal{F}(\tau) = \vert\tau^D\vert - \sigma_0$$
+with $\sigma_0$ as yield stress. If there is no plastic deformation - that is $\lambda=0$ - this yields $\vert\sigma^D\vert < \sigma_0$
+and otherwise if $\lambda > 0$ it follows that $\vert\sigma^D\vert = \sigma_0$. That means if the stress is smaller as the yield stress
+there are only elastic deformations. Therein the Index $D$ denotes the deviator part of the stress $\sigma$ which
+is dedined as
+$$\sigma^D = \sigma - \dfrac{1}{3}tr(\sigma).$$
+It describes the hydrostatic part of the stress tensor in contrast to the volumetric part. For metal the hydrostatic
+stress composes the main indicator for the plastic deformation.\\
+The second equation is called equilibrium condition with a force of areal density $f$ which we will neglect in our example.
+The boundary of $\Omega$ separates as follows $\Gamma=\Gamma_D\bigcup\Gamma_C$ and $\Gamma_D\bigcap\Gamma_C=\emptyset$.
+At the boundary $\Gamma_D$ we have zero Dirichlet conditions. $\Gamma_C$ denotes the potential contact boundary.\\
+The last two lines decribe the so-called Signorini contact conditions. If there is no contact the normal stress
+$$ \sigma_n = \sigma n\cdot n$$
+is zero with the outward normal $n$. If there is contact ($u_n = g$) the tangential stress $\sigma_t = \sigma\cdot n - \sigma_n n$
+vanishes, because we consider a frictionless situation and the normal stress is negative.
+
+\section{Derivation of the variational inequality}
+
+As a starting point we want to minimise an energy functional:
+$$E(\tau) := \dfrac{1}{2}\int\limits_{\Omega}\tau A \tau d\tau,\quad \tau\in \Pi W^{div}$$
+with
+$$W^{div}:=\lbrace \tau\in L^2(\Omega,\mathbb{R}^{dim\times\dim}_{sym}),div(\tau)\in L^2(\Omega,\mathbb{R}^{dim})\rbrace$$
+and
+$$\Pi \Sigma:=\lbrace \tau\in \Sigma, \mathcal{F}(\tau)\leq 0\rbrace$$
+as the set of admissible stresses which is defined
+by a continious, convex flow function $\mathcal{F}$.
+
+With the goal to derive the dual formulation of the minimisation problem, we define a lagrange function:
+$$L(\tau,\varphi) := E(\tau) + (\varphi, div(\tau)),\quad \lbrace\tau,\varphi\rbrace\in\Pi W^{div}\times U$$
+with $U := \lbrace u\in H^1(\Omega), u = g \text{ on } \Gamma_D,u_n\leq 0 \text{ on } \Gamma_C \rbrace$.\\
+By building the fr\'echet derivatives of $L$ for both components we obtain the dual formulation for the stationary case
+which is known as \textbf{Hencky-Type-Model}:\\
+Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times U$ with
+$$\left(A\sigma,\tau - \sigma\right) + \left(u, div(\tau) - div(\sigma)\right) \geq 0,\quad \forall \tau\in \Pi W^{div}$$
+$$-\left(div(\sigma),\varphi - u\right) \geq 0,\quad \forall \varphi\in U.$$
+By integrating by parts and multiplying the first inequality by $C=A^{-1}$ we achieve the primal-mixed version of our problem:\\
+Find a pair $\lbrace\sigma,u\rbrace\in \Pi W\times U$ with
+$$\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W$$
+$$\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in U.$$
+Therein $\varepsilon$ denotes the linearised deformation tensor with $\varepsilon(u) := \dfrac{1}{2}\left(\nabla u + \nabla u^T\right)$ for small deformations.\\
+Most materials - especially metals - have the property that they show some hardening effects during the forming process.
+There are different constitutive laws to describe those material behaviour. The most simple one is called linear isotropic hardening
+with the flow function $\mathcal{F}(\tau,\eta) = \vert\tau^D\vert - (\sigma_0 + \gamma\eta)$.
+It can be considered by establishing an additional term in our primal-mixed formulation:\\
+Find a pair $\lbrace(\sigma,\xi),u\rbrace\in \Pi (W\times L^2(\Omega,\mathbb{R}))\times U$ with
+$$\left(\sigma,\tau - \sigma\right) - \left(C\varepsilon(u), \tau - \sigma\right) + \gamma\left( \xi, \eta - \xi\right) \geq 0,\quad \forall (\tau,\eta)\in \Pi (W,L^2(\Omega,\mathbb{R}))$$
+$$\left(\sigma,\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in U,$$
+with the hardening parameter $\gamma > 0$.\\
+Now we want to derive a primal problem which only depends on the displacement $u$. For that purpose we
+set $\eta = \xi$ and eliminate the stress $\sigma$ by applying the projection theorem on\\
+$$\left(\sigma - C\varepsilon(u), \tau - \sigma\right) \geq 0,\quad \forall \tau\in \Pi W,$$
+which yields with the second inequality:\\
+Find the displacement $u\in U$ with
+$$\left(P_{\Pi}(C\varepsilon(u)),\varepsilon(\varphi) - \varepsilon(u)\right) \geq 0,\quad \forall \varphi\in U,$$
+with the projection:
+$$P_{\Pi}(\tau):=\begin{cases}
+ \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0 + \gamma\xi,\\
+ \hat\alpha\dfrac{\tau^D}{\vert\tau^D\vert} + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert > \sigma_0 + \gamma\xi,
+ \end{cases}$$
+with the radius
+$$\hat\alpha := \sigma_0 + \gamma\xi .$$
+With the relation $\xi = \vert\varepsilon(u) - A\sigma\vert$ it is possible to eliminate $\xi$ inside the projection $P_{\Pi}$:\\
+$$P_{\Pi}(\tau):=\begin{cases}
+ \tau, & \text{if }\vert\tau^D\vert \leq \sigma_0,\\
+ \alpha\dfrac{\tau^D}{\vert\tau^D\vert} + \dfrac{1}{3}tr(\tau), & \text{if }\vert\tau^D\vert > \sigma_0,
+ \end{cases}$$
+$$\alpha := \sigma_0 + \dfrac{\gamma}{2\mu+\gamma}\left(\vert\tau^D\vert - \sigma_0\right) ,$$
+with a further material parameter $\mu>0$ called shear modulus.\\
+So what we do is to calculate the stresses by using Hooke's law for linear elastic, isotropic materials
+$$\sigma = C \varepsilon(u) = 2\mu \varepsilon^D(u) + \kappa tr(\varepsilon(u))I = \left[2\mu\left(\mathbb{I} -\dfrac{1}{3} I\otimes I\right) + \kappa I\otimes I\right]\varepsilon(u)$$
+with the new material parameter $\kappa>0$ (bulk modulus). The variables $I$ and $\mathbb{I}$ denote the identity tensors of second and forth order.\\
+In the next step we test in a pointwise sense where the deviator part of the stress in a norm is bigger as the yield stress.
+If there are such points we project the deviator stress in those points back to the yield surface. Methods of this kind
+are called projections algorithm or radial-return-algorithm.\\
+Now we have a primal formulation of our elasto-plastic contact problem which only depends on the displacement $u$.
+It consists of a nonlinear variational inequality and has a unique solution as it shows the theorem of Lions and Stampaccia
+(A proof can be found in Rodrigues: Obstacle Problems in Mathematical Physics, North-Holland, Amsterdam, 1987).\\
+To handle the nonlinearity of the constitutive law we use a newton method and to deal with the contact we apply an
+active set method like in step-41. To be more concrete we combine both methods to an inexact semi smooth newton
+method - inexact since we use an iterative solver for the linearised problems in each newton step.
+
+\section{Linearisation of the constitutive law for the newton method}
+
+For the newton method we have to linearise the following semi-linearform
+$$a(\psi;\varphi) := \left(P_{\Pi}(C\varepsilon(\varphi)),\varepsilon(\varphi)\right).$$
+Becaus we have to find the solution $u$ in the convex set $U$, we have to apply an SQP-method (SQP: sequential quadratic
+programming). That means we have to solve a minimisation problem for a known $u^i$ in every SQP-step of the form
+\begin{eqnarray*}
+ & & a(u^{i};u^{i+1} - u^i) + \dfrac{1}{2}a'(u^i;u^{i+1} - u^i,u^{i+1} - u^i)\\
+ &=& a(u^i;u^{i+1}) - a(u^i;u^i) +\\
+ & & \dfrac{1}{2}\left( a'(u^i;u^{i+1},u^{i+1}) - 2a'(u^i;u^i,u^{i+1}) - a'(u^i;u^i,u^i)\right)\\
+ &\rightarrow& min,\quad u^{i+1}\in U.
+\end{eqnarray*}
+Neglecting the constant terms $ a(u^i;u^i)$ and $ a'(u^i;u^i,u^i)$ we obtain the following minimisation problem
+$$\dfrac{1}{2} a'(u^i;u^{i+1},u^{i+1}) - F(u^i)\rightarrow min,\quad u^{i+1}\in U$$
+with
+$$F(\varphi) := \left(a'(\varphi;\varphi,u^{i+1}) - a(\varphi;u^{i+1}) \right).$$
+In the case of our constitutive law the derivitive of the semi-linearform $a(.;.)$ at the point $u^i$ is
+
+$$a'(u^i;\psi,\varphi) =$$
+$$
+\begin{cases}
+\left(\left[2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes I\right) + \kappa I\otimes I\right]\varepsilon(\psi),\varepsilon(\varphi)\right), & \quad \vert\tau^D\vert \leq \sigma_0\\
+\left(\left[\dfrac{\alpha}{\vert\tau^D\vert}2\mu\left(\mathbb{I} - \dfrac{1}{3} I\otimes I - \dfrac{\tau^D\otimes\tau^D}{\vert\tau^D\vert}\right) + \kappa I\otimes I\right]\varepsilon(\psi),\varepsilon(\varphi) \right), & \quad \vert\tau^D\vert > \sigma_0
+\end{cases}
+$$
+with
+$$\tau^D := C\varepsilon^D(u^i).$$
+Again the first case is for elastic and the second for plastic deformation.
+
+\section{Formulation as a saddle point problem}
+
+On the line of step-41 we compose a saddle point problem out of the minimisation problem. Again we do so to gain a formulation
+that allows us to solve a linear system of equations finally.
+
+\section{Active Set methods to solve the saddle point problem}
+
+\section{The primal-dual active set algorithm combined with the inexact semi smooth newton method}
+
+The inexact newton method works as follows:
+\begin{itemize}
+ \item[(0)] Initialize $\mathcal{A}_k$ and $\mathcal{F}_k$, such that $\mathcal{S} = \mathcal{A}_k \cup \mathcal{F}_k$ and $\mathcal{A}_k \cap \mathcal{F}_k = \emptyset$ and set $k = 1$.
+ \item[(1)] Assembel the newton matrix $a'(U^k;\varphi_i,\varphi_j)$ and the right-hand-side $F(U^k)$.
+ \item[(2)] Find the primal-dual pair $(U^k,\Lambda^k)$ that satisfies
+ \begin{align*}
+ AU^k + B\Lambda^k & = F, &\\
+ \left[BU^k\right]_i & = G & & \forall i\in\mathcal{A}_k\\
+ \Lambda^k_i & = 0 & & \forall i\in\mathcal{F}_k.
+ \end{align*}
+% Note that $\mathcal{S}$ contains only dofs related to the boundary $\Gamma_C$. So in contrast to step-41 there are much more than $\vert \mathcal{S}\vert$ equations necessary to determine $U$ and $\Lambda$.
+ \item[(3)] Define the new active and inactive sets by
+ $$\mathcal{A}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c\left(\left[BU^k\right]_i - G_i\left) < 0\rbrace,$$
+ $$\mathcal{F}_{k+1}:=\lbrace i\in\mathcal{S}:\Lambda^k_i + c\left(\left[BU^k\right]_i - G_i\left) \geq 0\rbrace.$$
+ \item[(4)] If $\mathcal{A}_{k+1} = \mathcal{A}_k$ and $\vert F(U^{k+1}\vert < \delta$ then stop, else set $k=k+1$ and go to step (1).
+\end{itemize}
+
+\section{Implementation}
+
+\end{document}
\ No newline at end of file
}
}
}
+
+/**
+ * Laplacian residual operator in weak form
+ *
+ * \f[
+ * -(\nabla u, \nabla v)
+ * \f]
+ */
+ template <int dim>
+ inline void
+ cell_residual (
+ Vector<double>& result,
+ const FEValuesBase<dim>& fe,
+ const std::vector<Tensor<1,dim> >& input,
+ double factor = 1.)
+ {
+ const unsigned int nq = fe.n_quadrature_points;
+ const unsigned int n_dofs = fe.dofs_per_cell;
+ Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq));
+ Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
+
+ for (unsigned k=0;k<nq;++k)
+ {
+ const double dx = factor * fe.JxW(k);
+ for (unsigned i=0;i<n_dofs;++i)
+ result(i) += dx * (input[k] * fe.shape_grad(i,k));
+ }
+ }
+
+
+/**
+ * Vector-valued Laplacian residual operator in weak form
+ *
+ * \f[
+ * -(\nabla u, \nabla v)
+ * \f]
+ */
+ template <int dim>
+ inline void
+ cell_residual (
+ Vector<double>& result,
+ const FEValuesBase<dim>& fe,
+ const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& input,
+ double factor = 1.)
+ {
+ const unsigned int nq = fe.n_quadrature_points;
+ const unsigned int n_dofs = fe.dofs_per_cell;
+ const unsigned int n_comp = fe.get_fe().n_components();
+
+ AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
+ Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs));
+
+ for (unsigned k=0;k<nq;++k)
+ {
+ const double dx = factor * fe.JxW(k);
+ for (unsigned i=0;i<n_dofs;++i)
+ for (unsigned int d=0;d<n_comp;++d)
+ {
+
+ result(i) += dx * (input[d][k] * fe.shape_grad_component(i,k,d));
+ }
+ }
+ }
+
+
/**
* Weak boundary condition of Nitsche type for the Laplacian, namely on the face <i>F</i> the matrix
* @f[
}
/**
- * Weak boundary condition for the Laplace operator by Nitsche, vector
- * valued version, namely on the face <i>F</i>
- * the vector
+ * Weak boundary condition for the Laplace operator by Nitsche, scalar
+ * version, namely on the face <i>F</i> the vector
* @f[
* \int_F \Bigl(\gamma (u-g) v - \partial_n u v - (u-g) \partial_n v\Bigr)\;ds.
* @f]
void nitsche_residual (
Vector<double>& result,
const FEValuesBase<dim>& fe,
- const VectorSlice<const std::vector<std::vector<double> > >& input,
- const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Dinput,
- const VectorSlice<const std::vector<std::vector<double> > >& data,
+ const std::vector<double>& input,
+ const std::vector<Tensor<1,dim> >& Dinput,
+ const std::vector<double>& data,
double penalty,
double factor = 1.)
{
const unsigned int n_dofs = fe.dofs_per_cell;
-
- const unsigned int n_comp = fe.get_fe().n_components();
- AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
- AssertVectorVectorDimension(Dinput, n_comp, fe.n_quadrature_points);
- AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);
+ AssertDimension(input.size(), fe.n_quadrature_points);
+ AssertDimension(Dinput.size(), fe.n_quadrature_points);
+ AssertDimension(data.size(), fe.n_quadrature_points);
for (unsigned k=0;k<fe.n_quadrature_points;++k)
{
const double dx = factor * fe.JxW(k);
const Point<dim>& n = fe.normal_vector(k);
for (unsigned i=0;i<n_dofs;++i)
- for (unsigned int d=0;d<n_comp;++d)
- {
- const double dnv = fe.shape_grad_component(i,k,d) * n;
- const double dnu = Dinput[d][k] * n;
- const double v= fe.shape_value_component(i,k,d);
- const double u= input[d][k];
- const double g= data[d][k];
+ {
+ const double dnv = fe.shape_grad(i,k) * n;
+ const double dnu = Dinput[k] * n;
+ const double v= fe.shape_value(i,k);
+ const double u= input[k];
+ const double g= data[k];
- result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
- }
+ result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
+ }
}
}
/**
- * Weak boundary condition for the Laplace operator by Nitsche, scalar
- * version, namely on the face <i>F</i> the vector
+ * Weak boundary condition for the Laplace operator by Nitsche, vector
+ * valued version, namely on the face <i>F</i>
+ * the vector
* @f[
* \int_F \Bigl(\gamma (u-g) v - \partial_n u v - (u-g) \partial_n v\Bigr)\;ds.
* @f]
* <tt>Dinput</tt>, respectively. <i>g</i> is the inhomogeneous
* boundary value in the argument <tt>data</tt>. $\gamma$ is the usual
* penalty parameter.
+ *
+ * @ingroup Integrators
+ * @author Guido Kanschat
+ * @date 2008, 2009, 2010
*/
template <int dim>
void nitsche_residual (
Vector<double>& result,
const FEValuesBase<dim>& fe,
- const std::vector<double>& input,
- const std::vector<Tensor<1,dim> >& Dinput,
- const std::vector<double>& data,
+ const VectorSlice<const std::vector<std::vector<double> > >& input,
+ const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Dinput,
+ const VectorSlice<const std::vector<std::vector<double> > >& data,
double penalty,
double factor = 1.)
{
const unsigned int n_dofs = fe.dofs_per_cell;
- AssertDimension(input.size(), fe.n_quadrature_points);
- AssertDimension(Dinput.size(), fe.n_quadrature_points);
- AssertDimension(data.size(), fe.n_quadrature_points);
+ const unsigned int n_comp = fe.get_fe().n_components();
+ AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points);
+ AssertVectorVectorDimension(Dinput, n_comp, fe.n_quadrature_points);
+ AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points);
for (unsigned k=0;k<fe.n_quadrature_points;++k)
{
const double dx = factor * fe.JxW(k);
const Point<dim>& n = fe.normal_vector(k);
for (unsigned i=0;i<n_dofs;++i)
- {
- const double dnv = fe.shape_grad(i,k) * n;
- const double dnu = Dinput[k] * n;
- const double v= fe.shape_value(i,k);
- const double u= input[k];
- const double g= data[k];
+ for (unsigned int d=0;d<n_comp;++d)
+ {
+ const double dnv = fe.shape_grad_component(i,k,d) * n;
+ const double dnu = Dinput[d][k] * n;
+ const double v= fe.shape_value_component(i,k,d);
+ const double u= input[d][k];
+ const double g= data[d][k];
- result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
- }
+ result(i) += dx*(2.*penalty*(u-g)*v - dnv*(u-g) - dnu*v);
+ }
}
}
}
}
}
+
+/**
+ * Residual term for the symmetric interior penalty method.
+ *
+ * @ingroup Integrators
+ * @author Guido Kanschat
+ * @date 2012
+ */
+ template<int dim>
+ void
+ ip_residual(
+ Vector<double>& result1,
+ Vector<double>& result2,
+ const FEValuesBase<dim>& fe1,
+ const FEValuesBase<dim>& fe2,
+ const std::vector<double>& input1,
+ const std::vector<Tensor<1,dim> >& Dinput1,
+ const std::vector<double>& input2,
+ const std::vector<Tensor<1,dim> >& Dinput2,
+ double pen,
+ double int_factor = 1.,
+ double ext_factor = -1.)
+{
+ Assert(fe1.get_fe().n_components() == 1,
+ ExcDimensionMismatch(fe1.get_fe().n_components(), 1));
+ Assert(fe2.get_fe().n_components() == 1,
+ ExcDimensionMismatch(fe2.get_fe().n_components(), 1));
+
+ const double nui = int_factor;
+ const double nue = (ext_factor < 0) ? int_factor : ext_factor;
+ const double penalty = .5 * pen * (nui + nue);
+
+ const unsigned int n_dofs = fe1.dofs_per_cell;
+
+ for (unsigned k=0;k<fe1.n_quadrature_points;++k)
+ {
+ const double dx = fe1.JxW(k);
+ const Point<dim>& n = fe1.normal_vector(k);
+
+ for (unsigned i=0;i<n_dofs;++i)
+ {
+ const double vi = fe1.shape_value(i,k);
+ const Tensor<1,dim>& Dvi = fe1.shape_grad(i,k);
+ const double dnvi = Dvi * n;
+ const double ve = fe2.shape_value(i,k);
+ const Tensor<1,dim>& Dve = fe2.shape_grad(i,k);
+ const double dnve = Dve * n;
+
+ const double ui = input1[k];
+ const Tensor<1,dim>& Dui = Dinput1[k];
+ const double dnui = Dui * n;
+ const double ue = input2[k];
+ const Tensor<1,dim>& Due = Dinput2[k];
+ const double dnue = Due * n;
+
+ result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
+ result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
+ result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
+ result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
+ }
+ }
+}
+
+
+/**
+ * Vector-valued residual term for the symmetric interior penalty method.
+ *
+ * @ingroup Integrators
+ * @author Guido Kanschat
+ * @date 2012
+ */
+ template<int dim>
+ void
+ ip_residual(
+ Vector<double>& result1,
+ Vector<double>& result2,
+ const FEValuesBase<dim>& fe1,
+ const FEValuesBase<dim>& fe2,
+ const VectorSlice<const std::vector<std::vector<double> > >& input1,
+ const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Dinput1,
+ const VectorSlice<const std::vector<std::vector<double> > >& input2,
+ const VectorSlice<const std::vector<std::vector<Tensor<1,dim> > > >& Dinput2,
+ double pen,
+ double int_factor = 1.,
+ double ext_factor = -1.)
+ {
+ const unsigned int n_comp = fe1.get_fe().n_components();
+ const unsigned int n1 = fe1.dofs_per_cell;
+
+ AssertVectorVectorDimension(input1, n_comp, fe1.n_quadrature_points);
+ AssertVectorVectorDimension(Dinput1, n_comp, fe1.n_quadrature_points);
+ AssertVectorVectorDimension(input2, n_comp, fe2.n_quadrature_points);
+ AssertVectorVectorDimension(Dinput2, n_comp, fe2.n_quadrature_points);
+
+ const double nui = int_factor;
+ const double nue = (ext_factor < 0) ? int_factor : ext_factor;
+ const double penalty = .5 * pen * (nui + nue);
+
+
+ for (unsigned k=0;k<fe1.n_quadrature_points;++k)
+ {
+ const double dx = fe1.JxW(k);
+ const Point<dim>& n = fe1.normal_vector(k);
+
+ for (unsigned i=0;i<n1;++i)
+ for (unsigned int d=0;d<n_comp;++d)
+ {
+ const double vi = fe1.shape_value_component(i,k,d);
+ const Tensor<1,dim>& Dvi = fe1.shape_grad_component(i,k,d);
+ const double dnvi = Dvi * n;
+ const double ve = fe2.shape_value_component(i,k,d);
+ const Tensor<1,dim>& Dve = fe2.shape_grad_component(i,k,d);
+ const double dnve = Dve * n;
+
+ const double ui = input1[d][k];
+ const Tensor<1,dim>& Dui = Dinput1[d][k];
+ const double dnui = Dui * n;
+ const double ue = input2[d][k];
+ const Tensor<1,dim>& Due = Dinput2[d][k];
+ const double dnue = Due * n;
+
+ result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi);
+ result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue);
+ result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve);
+ result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve);
+ }
+ }
+ }
+
+
/**
* Auxiliary function computing the penalty parameter for interior
* <i>h<sub>i</sub></i> is the length of <i>Z<sub>i</sub></i>
* orthogonal to the current face.
*
+ * @ingroup Integrators
* @author Guido Kanschat
* @date 2010
*/
void Tvmult_add (VectorBase &dst,
const VectorBase &src) const;
- /**
- * Return the square of the norm
- * of the vector $v$ with respect
- * to the norm induced by this
- * matrix,
- * i.e. $\left(v,Mv\right)$. This
- * is useful, e.g. in the finite
- * element context, where the
- * $L_2$ norm of a function
- * equals the matrix norm with
- * respect to the mass matrix of
- * the vector representing the
- * nodal values of the finite
- * element function.
- *
- * Obviously, the matrix needs to
- * be quadratic for this operation.
- *
- * The implementation of this function
- * is not as efficient as the one in
- * the @p MatrixBase class used in
- * deal.II (i.e. the original one, not
- * the PETSc wrapper class) since PETSc
- * doesn't support this operation and
- * needs a temporary vector.
- *
- * Note that if the current object
- * represents a parallel distributed
- * matrix (of type
- * PETScWrappers::MPI::SparseMatrix),
- * then the given vector has to be
- * a distributed vector as
- * well. Conversely, if the matrix is
- * not distributed, then neither
- * may the vector be.
- */
- PetscScalar matrix_norm_square (const VectorBase &v) const;
-
- /**
- * Compute the matrix scalar
- * product $\left(u,Mv\right)$.
- *
- * The implementation of this function
- * is not as efficient as the one in
- * the @p MatrixBase class used in
- * deal.II (i.e. the original one, not
- * the PETSc wrapper class) since PETSc
- * doesn't support this operation and
- * needs a temporary vector.
- *
- * Note that if the current object
- * represents a parallel distributed
- * matrix (of type
- * PETScWrappers::MPI::SparseMatrix),
- * then both vectors have to be
- * distributed vectors as
- * well. Conversely, if the matrix is
- * not distributed, then neither of the
- * vectors may be.
- */
- PetscScalar matrix_scalar_product (const VectorBase &u,
- const VectorBase &v) const;
/**
* Compute the residual of an
# include <deal.II/lac/exceptions.h>
# include <deal.II/lac/petsc_matrix_base.h>
+# include <deal.II/lac/petsc_parallel_vector.h>
# include <vector>
DEAL_II_NAMESPACE_OPEN
namespace MPI
{
+
+
/**
* Implementation of a parallel sparse matrix class based on PETSC, with rows
* of the matrix distributed across an MPI network. All the functionality is
class SparseMatrix : public MatrixBase
{
public:
+
/**
* A structure that describes some of
* the traits of this class in terms
<< "The number of local rows " << arg1
<< " must be larger than the total number of rows " << arg2);
//@}
+
+ /**
+ * Return the square of the norm
+ * of the vector $v$ with respect
+ * to the norm induced by this
+ * matrix,
+ * i.e. $\left(v,Mv\right)$. This
+ * is useful, e.g. in the finite
+ * element context, where the
+ * $L_2$ norm of a function
+ * equals the matrix norm with
+ * respect to the mass matrix of
+ * the vector representing the
+ * nodal values of the finite
+ * element function.
+ *
+ * Obviously, the matrix needs to
+ * be quadratic for this operation.
+ *
+ * The implementation of this function
+ * is not as efficient as the one in
+ * the @p MatrixBase class used in
+ * deal.II (i.e. the original one, not
+ * the PETSc wrapper class) since PETSc
+ * doesn't support this operation and
+ * needs a temporary vector.
+ */
+ PetscScalar matrix_norm_square (const Vector &v) const;
+
+ /**
+ * Compute the matrix scalar
+ * product $\left(u,Mv\right)$.
+ *
+ * The implementation of this function
+ * is not as efficient as the one in
+ * the @p MatrixBase class used in
+ * deal.II (i.e. the original one, not
+ * the PETSc wrapper class) since PETSc
+ * doesn't support this operation and
+ * needs a temporary vector.
+ */
+ PetscScalar matrix_scalar_product (const Vector &u,
+ const Vector &v) const;
+
private:
/**
class SparseMatrix : public MatrixBase
{
public:
+
/**
* A structure that describes some of
* the traits of this class in terms of
*/
virtual const MPI_Comm & get_mpi_communicator () const;
+ /**
+ * Return the square of the norm
+ * of the vector $v$ with respect
+ * to the norm induced by this
+ * matrix,
+ * i.e. $\left(v,Mv\right)$. This
+ * is useful, e.g. in the finite
+ * element context, where the
+ * $L_2$ norm of a function
+ * equals the matrix norm with
+ * respect to the mass matrix of
+ * the vector representing the
+ * nodal values of the finite
+ * element function.
+ *
+ * Obviously, the matrix needs to
+ * be quadratic for this operation.
+ *
+ * The implementation of this function
+ * is not as efficient as the one in
+ * the @p MatrixBase class used in
+ * deal.II (i.e. the original one, not
+ * the PETSc wrapper class) since PETSc
+ * doesn't support this operation and
+ * needs a temporary vector.
+ */
+ PetscScalar matrix_norm_square (const VectorBase &v) const;
+
+ /**
+ * Compute the matrix scalar
+ * product $\left(u,Mv\right)$.
+ *
+ * The implementation of this function
+ * is not as efficient as the one in
+ * the @p MatrixBase class used in
+ * deal.II (i.e. the original one, not
+ * the PETSc wrapper class) since PETSc
+ * doesn't support this operation and
+ * needs a temporary vector.
+ */
+ PetscScalar matrix_scalar_product (const VectorBase &u,
+ const VectorBase &v) const;
+
private:
/**
const MPI_Comm &communicator = MPI_COMM_WORLD);
/**
- * Reinit functionality. This
- * function destroys the old vector
- * content and generates a new one
- * based on the input partitioning.
+ * Reinit functionality. This function
+ * destroys the old vector content and
+ * generates a new one based on the
+ * input partitioning. The flag
+ * <tt>fast</tt> determines whether the
+ * vector should be filled with zero
+ * (false) or left untouched (true).
*/
void reinit (const IndexSet ¶llel_partitioning,
const MPI_Comm &communicator = MPI_COMM_WORLD,
/**
* Initialization with an
- * IndexSet. Similar to the call in
- * the other class MPI::Vector, with
- * the difference that now a copy on
- * all processes is generated. This
+ * IndexSet. Similar to the call in the
+ * other class MPI::Vector, with the
+ * difference that now a copy on all
+ * processes is generated. This
* initialization function is
- * appropriate in case the data in
- * the localized vector should be
- * imported from a distributed vector
- * that has been initialized with the
- * same communicator. The variable
- * <tt>fast</tt> determines whether
- * the vector should be filled with
- * zero or left untouched.
+ * appropriate in case the data in the
+ * localized vector should be imported
+ * from a distributed vector that has
+ * been initialized with the same
+ * communicator. The variable
+ * <tt>fast</tt> determines whether the
+ * vector should be filled with zero
+ * (false) or left untouched (true).
*/
void reinit (const IndexSet &input_map,
const MPI_Comm &communicator = MPI_COMM_WORLD,
VectorBase &
VectorBase::operator = (const TrilinosScalar s)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (numbers::is_finite(s), ExcNumberNotFinite());
VectorBase::set (const std::vector<unsigned int> &indices,
const std::vector<TrilinosScalar> &values)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
Assert (indices.size() == values.size(),
ExcDimensionMismatch(indices.size(),values.size()));
VectorBase::set (const std::vector<unsigned int> &indices,
const ::dealii::Vector<TrilinosScalar> &values)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
+
Assert (indices.size() == values.size(),
ExcDimensionMismatch(indices.size(),values.size()));
void
VectorBase::add (const TrilinosScalar s)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (numbers::is_finite(s), ExcNumberNotFinite());
unsigned int n_local = local_size();
VectorBase::add (const TrilinosScalar a,
const VectorBase &v)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == v.local_size(),
ExcDimensionMismatch(local_size(), v.local_size()));
const TrilinosScalar b,
const VectorBase &w)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == v.local_size(),
ExcDimensionMismatch(local_size(), v.local_size()));
Assert (local_size() == w.local_size(),
VectorBase::sadd (const TrilinosScalar s,
const VectorBase &v)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == v.local_size(),
ExcDimensionMismatch(local_size(), v.local_size()));
const TrilinosScalar a,
const VectorBase &v)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == v.local_size(),
ExcDimensionMismatch(local_size(), v.local_size()));
const TrilinosScalar b,
const VectorBase &w)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == v.local_size(),
ExcDimensionMismatch(local_size(), v.local_size()));
Assert (local_size() == w.local_size(),
const TrilinosScalar c,
const VectorBase &x)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == v.local_size(),
ExcDimensionMismatch(local_size(), v.local_size()));
Assert (local_size() == w.local_size(),
void
VectorBase::scale (const VectorBase &factors)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (local_size() == factors.local_size(),
ExcDimensionMismatch(local_size(), factors.local_size()));
VectorBase::equ (const TrilinosScalar a,
const VectorBase &v)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (numbers::is_finite(a), ExcNumberNotFinite());
// If we don't have the same map, copy.
const TrilinosScalar b,
const VectorBase &w)
{
+ // if we have ghost values, do not allow
+ // writing to this vector at all.
+ Assert (!has_ghost_elements(), ExcGhostsPresent());
Assert (v.local_size() == w.local_size(),
ExcDimensionMismatch (v.local_size(), w.local_size()));
* course at the expense of getting only an approximate result.)
* </ol>
*
- * When recording a new mnemonic name, the user must supply a std::vector<unsigned int>
- * components that lists the @ref GlossComponent "(vector) components" to be extracted from the
- * given input. If the user simply wants to extract all the components, or the first
- * n components then default parameters and overloaded functions mean that
- * the std::vector need not be explicitly supplied to the @p add_field_name method.
- * If the @p evaluate_field with a @p DataPostprocessor object is used, the components
- * are interpreted as the component of the @p DataPostprocessor return vector. This
- * vector can be larger than the FE space.
+ * When recording a new mnemonic name, the user must supply a
+ * std::vector@<bool@> component_mask (see @ref GlossComponentMask "this glossary entry")
+ * to indicate the @ref GlossComponent "(vector) components"
+ * to be extracted from the given input. If the user simply wants to extract
+ * all the components, the mask need not be explicitly supplied to the @p
+ * add_field_name method and the default value of the parameter is sufficient.
+ * If the @p evaluate_field with a @p DataPostprocessor object is used, the
+ * component_mask is interpreted as the mask of the @p DataPostprocessor
+ * return vector. The size of this mask can be different to that of the FE
+ * space, but must be provided when the @p add_field_name method is
+ * called. One variant of the @p add_field_name method allows an unsigned int
+ * input to construct a suitable mask, if all values from the @p
+ * DataPostprocessor are desired.
*
* The class automatically generates names for the data stored based on the mnemonics
* supplied. The methods @p add_component_names and @p add_independent_names allow
/**
* Put another mnemonic string (and hence
* @p VECTOR) into the class. This method
- * adds storage space for components.size()
- * variables.
+ * adds storage space for variables equal
+ * to the number of true values in
+ * component_mask.
* This also adds extra entries for points
* that are already in the class, so
* @p add_field_name and @p add_points can
* be called in any order.
*/
void add_field_name(const std::string &vector_name,
- const std::vector <unsigned int> &components = std::vector <unsigned int>());
+ const std::vector <bool> &component_mask = std::vector <bool>());
/**
* Put another mnemonic string (and hence
* @p VECTOR) into the class. This method
- * adds storage space for components.size()
+ * adds storage space for n_components
* variables.
* This also adds extra entries for points
* that are already in the class, so
void add_field_name(const std::string &vector_name,
const unsigned int n_components);
- /**
- * Provide optional names for each component
- * of a field. These names will be used
- * instead of names generated from the
- * field name, if supplied.
- */
+ /**
+ * Provide optional names for each component
+ * of a field. These names will be used
+ * instead of names generated from the
+ * field name, if supplied.
+ */
void add_component_names(const std::string &vector_name,
- const std::vector <std::string> &component_names);
-
- /**
- * Provide optional names for the
- * independent values. These names will
- * be used instead of "Indep_...", if
- * supplied.
- */
+ const std::vector <std::string> &component_names);
+
+ /**
+ * Provide optional names for the
+ * independent values. These names will
+ * be used instead of "Indep_...", if
+ * supplied.
+ */
void add_independent_names(const std::vector <std::string> &independent_names);
* Extract values at the stored points
* from the VECTOR supplied and add them
* to the new dataset in vector_name.
- * The list of components supplied when the
+ * The component mask supplied when the
* field was added is used to select
* components to extract.
* If a @p DoFHandler is used, one (and only
* Compute values using a @p DataPostprocessor object
* with the @p VECTOR supplied and add them
* to the new dataset in vector_name.
- * The list of components supplied when the
+ * The component_mask supplied when the
* field was added is used to select
* components to extract from the @p DataPostprocessor
* return vector.
* data. Therefore, if only this method is used,
* the class is fully compatible with
* adaptive refinement.
- * The list of components supplied when the
+ * The component_mask supplied when the
* field was added is used to select
* components to extract. If
* a @p DoFHandler is used, one (and only
* Check the internal data sizes to test
* for a loss of data sync. This is often
* used in @p Assert statements with the
- * @p ExcDataLostSync exception. If @p
- * strict is @p false this method returns
+ * @p ExcDataLostSync exception. If @p strict
+ * is @p false this method returns
* @p true if all sizes are within 1 of
* each other (needed to allow data to be
* added), with @p strict = @p true they
DeclException0(ExcDoFHandlerRequired);
/**
- * The triangulation indicated that mesh
+ * The triangulation indicated that mesh
* has been refined in some way. This suggests
* that the internal dof indices stored
* are no longer meaningful.
*/
std::map <std::string, std::vector <std::vector <double> > > data_store;
+ /**
+ * Saves a component mask
+ * for each mnemonic.
+ */
+ std::map <std::string, std::vector<bool> > component_mask;
- /**
- * Saves a vector listing components
- * associated with a mnemonic.
- */
- std::map <std::string, std::vector<unsigned int> > field_components;
/**
* Saves a vector listing component
unsigned int n_indep;
- /**
- * A function that will be triggered
- * through signals whenever the
- * triangulation is modified.
- *
- * It is currently used to check
- * if the triangulation has changed,
- * invalidating precomputed values.
- */
+ /**
+ * A function that will be triggered
+ * through signals whenever the
+ * triangulation is modified.
+ *
+ * It is currently used to check
+ * if the triangulation has changed,
+ * invalidating precomputed values.
+ */
void tria_change_listener ();
};
}
-
- PetscScalar
- MatrixBase::matrix_norm_square (const VectorBase &v) const
- {
- Vector tmp(v.size());
- vmult (tmp, v);
- return tmp*v;
- }
-
-
-
- PetscScalar
- MatrixBase::matrix_scalar_product (const VectorBase &u,
- const VectorBase &v) const
- {
- Vector tmp(v.size());
- vmult (tmp, v);
- return u*tmp;
- }
-
-
-
PetscScalar
MatrixBase::residual (VectorBase &dst,
const VectorBase &x,
}
}
-
-
// explicit instantiations
//
template
const std::vector<unsigned int> &,
const unsigned int ,
const bool);
+
+
+ PetscScalar
+ SparseMatrix::matrix_norm_square (const Vector &v) const
+ {
+ Vector tmp (v);
+ vmult (tmp, v);
+ return tmp*v;
+ }
+
+ PetscScalar
+ SparseMatrix::matrix_scalar_product (const Vector &u,
+ const Vector &v) const
+ {
+ Vector tmp (v);
+ vmult (tmp, v);
+ return u*tmp;
+ }
+
}
}
namespace PETScWrappers
{
+
SparseMatrix::SparseMatrix ()
{
const int m=0, n=0, n_nonzero_per_row=0;
template void
SparseMatrix::do_reinit (const CompressedSimpleSparsityPattern &,
const bool);
+
+ PetscScalar
+ SparseMatrix::matrix_norm_square (const VectorBase &v) const
+ {
+ Vector tmp (v.size());
+ vmult (tmp, v);
+ return tmp*v;
+ }
+
+ PetscScalar
+ SparseMatrix::matrix_scalar_product (const VectorBase &u,
+ const VectorBase &v) const
+ {
+ Vector tmp (v.size());
+ vmult (tmp, v);
+ return u*tmp;
+ }
}
namespace Algorithms
{
-#include "operator.inst"
+#include "dof_output_operator.inst"
}
DEAL_II_NAMESPACE_CLOSE
#include <deal.II/numerics/point_value_history.h>
+#include <algorithm>
+
DEAL_II_NAMESPACE_OPEN
independent_values = point_value_history.independent_values;
indep_names = point_value_history.indep_names;
data_store = point_value_history.data_store;
- field_components = point_value_history.field_components;
+ component_mask = point_value_history.component_mask;
component_names_map = point_value_history.component_names_map;
point_geometry_data = point_value_history.point_geometry_data;
independent_values = point_value_history.independent_values;
indep_names = point_value_history.indep_names;
data_store = point_value_history.data_store;
- field_components = point_value_history.field_components;
+ component_mask = point_value_history.component_mask;
component_names_map = point_value_history.component_names_map;
point_geometry_data = point_value_history.point_geometry_data;
// check each cell to find a suitable
// support points
+ // GridTools::find_active_cell_around_point
+ // is an alternative. That method is not
+ // used here mostly because of the history
+ // of the class. The algorithm used in
+ // add_points below may be slightly more
+ // efficient than find_active_cell_around_point
+ // because it operates on a set of points.
+
for (; cell != endc; cell++)
{
fe_values.reinit (cell);
// will still lie within one cell distance.
// calling
// GridTools::find_active_cell_around_point
- // to obtain a cell to search may be an
+ // to obtain a cell to search is an
// option for these methods, but currently
// the GridTools method does not cater for
// a vector of points, and does not seem to
{
// add an extra row to each vector
// entry
- unsigned int n_components = (field_components.find (data_store_begin->first))->second.size();
- for (unsigned int component = 0; component < n_components; component++)
+ std::vector <bool> current_mask = (component_mask.find (data_store_begin->first))->second;
+ unsigned int n_stored = std::count(current_mask.begin(), current_mask.end(), true);
+ for (unsigned int component = 0; component < n_stored; component++)
{
data_store_begin->second.push_back (std::vector<double> (0));
}
// check each cell to find suitable support
// points
+ // GridTools::find_active_cell_around_point
+ // is an alternative. That method is not
+ // used here mostly because of the history
+ // of the class. The algorithm used here
+ // may be slightly more
+ // efficient than find_active_cell_around_point
+ // because it operates on a set of points.
for (; cell != endc; cell++)
{
fe_values.reinit (cell);
{
// add an extra row to each vector
// entry
- unsigned int n_components = (field_components.find (data_store_begin->first))->second.size();
- for (unsigned int component = 0; component < n_components; component++)
+ std::vector <bool> current_mask = (component_mask.find (data_store_begin->first))->second;
+ unsigned int n_stored = std::count(current_mask.begin(), current_mask.end(), true);
+ for (unsigned int component = 0; component < n_stored; component++)
{
data_store_begin->second.push_back (std::vector<double> (0));
}
template <int dim>
void PointValueHistory<dim>
-::add_field_name (const std::string &vector_name, const std::vector <unsigned int> &components)
+::add_field_name (const std::string &vector_name, const std::vector <bool> &mask)
{
// can't be closed to add additional points
// or vectors
AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ());
- // Make a field_components list, if not supplied
- std::vector <unsigned int> temp_components = components;
- if (temp_components.size() == 0)
+ // Make a component_mask, if not supplied
+ std::vector <bool> temp_mask = mask;
+ if (temp_mask.size() == 0)
{
- // Generate component list
- // 0,..,fe.n_components -1
- for (unsigned int component = 0; component < dof_handler->get_fe ().n_components (); component++)
- temp_components.push_back (component);
+ temp_mask = std::vector<bool> (dof_handler->get_fe().n_components(), true);
}
- field_components.insert (std::pair <std::string, std::vector<unsigned int> > (vector_name, temp_components));
+ component_mask.insert (std::pair <std::string, std::vector<bool> > (vector_name, temp_mask));
// insert an empty vector of strings
// to ensure each field has an entry
// point_geometry_data.size() long
std::pair<std::string, std::vector <std::vector <double> > > pair_data;
pair_data.first = vector_name;
- int n_datastreams = point_geometry_data.size () * temp_components.size(); // each point has field_components sub parts
+ unsigned int n_stored = std::count(temp_mask.begin(), temp_mask.end(), true);
+ int n_datastreams = point_geometry_data.size () * n_stored; // each point has n_stored sub parts
std::vector < std::vector <double> > vector_size (n_datastreams,
std::vector <double> (0));
pair_data.second = vector_size;
void PointValueHistory<dim>
::add_field_name(const std::string &vector_name, const unsigned int n_components)
{
- std::vector <unsigned int> temp_components;
- // Generate component list
- // 0,..,n_components -1
- for (unsigned int component = 0; component < n_components; component++)
- temp_components.push_back (component);
-
- add_field_name (vector_name, temp_components);
+ std::vector <bool> temp_mask (n_components, true);
+ add_field_name (vector_name, temp_mask);
}
typename std::map <std::string, std::vector <std::string> >::iterator names = component_names_map.find(vector_name);
Assert (names != component_names_map.end(), ExcMessage("vector_name not in class"));
- typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
- Assert (components != field_components.end(), ExcMessage("vector_name not in class"));
- Assert (component_names.size() == components->second.size(), ExcDimensionMismatch (component_names.size(), components->second.size()));
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(vector_name);
+ Assert (mask != component_mask.end(), ExcMessage("vector_name not in class"));
+ unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true);
+ Assert (component_names.size() == n_stored, ExcDimensionMismatch (component_names.size(), n_stored));
names->second = component_names;
}
// to check vector_name is in the map.
typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(vector_name);
Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
- // Repeat for field_components
- typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
- Assert (components != field_components.end(), ExcMessage("vector_name not in class"));
+ // Repeat for component_mask
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(vector_name);
+ Assert (mask != component_mask.end(), ExcMessage("vector_name not in class"));
- Assert (components->second.size () <= dof_handler->get_fe ().n_components (), ExcIndexRange (components->second.size (), 0, dof_handler->get_fe ().n_components () + 1)); // 0 components legal as request to generate list internally!
- for (unsigned int component = 0; component < components->second.size (); component++)
- {
- Assert (components->second [component] < dof_handler->get_fe ().n_components (), ExcIndexRange (components->second[component], 0, dof_handler->get_fe ().n_components ()));
- } // hopefully the compiler will optimise
+ Assert (mask->second.size () == dof_handler->get_fe ().n_components (), ExcDimensionMismatch (mask->second.size (), dof_handler->get_fe ().n_components ()));
+ unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true);
typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
for (unsigned int data_store_index = 0; point != point_geometry_data.end (); point++, data_store_index++)
{
- // Look up the component to add
- // in field_components, and
+ // Look up the components to add
+ // in the component_mask, and
// access the data associated with
- // that component
- for (unsigned int comp = 0; comp < components->second.size(); comp++)
+ // those components
+
+ for (unsigned int store_index = 0, comp = 0; comp < mask->second.size(); comp++)
{
- unsigned int solution_index = point->solution_indices[components->second[comp]];
- data_store_field->second[data_store_index * components->second.size() + comp].push_back (solution (solution_index));
+ if (mask->second[comp])
+ {
+ unsigned int solution_index = point->solution_indices[comp];
+ data_store_field->second[data_store_index * n_stored + store_index].push_back (solution (solution_index));
+ store_index++;
+ }
}
}
}
{
typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(*name);
Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
- // Repeat for field_components
- typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(*name);
- Assert (components != field_components.end(), ExcMessage("vector_name not in class"));
-
- // Push back components.size() computed
- // quantities. Each component requested
- // must be less than the postprocessor
- // n_output_variables.
- for (unsigned int comp = 0; comp < components->second.size(); comp++)
+ // Repeat for component_mask
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(*name);
+ Assert (mask != component_mask.end(), ExcMessage("vector_name not in class"));
+
+ Assert (mask->second.size () == n_output_variables, ExcDimensionMismatch (mask->second.size (),n_output_variables));
+ unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true);
+
+ // Push back computed quantities according
+ // to the component_mask.
+ for (unsigned int store_index = 0, comp = 0; comp < mask->second.size(); comp++)
{
- Assert (components->second[comp] < n_output_variables, ExcIndexRange (components->second[comp], 0, n_output_variables));
- data_store_field->second[data_store_index * components->second.size() + comp].push_back (computed_quantities[0](components->second[comp]));
+ if (mask->second[comp])
+ {
+ data_store_field->second[data_store_index * n_stored + store_index].push_back (computed_quantities[0](comp));
+ store_index++;
+ }
}
}
} // end of loop over points
AssertThrow (have_dof_handler, ExcDoFHandlerRequired ());
if (n_indep != 0) // hopefully this will get optimized, can't test independent_values[0] unless n_indep > 0
- {
- Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
- }
+ {
+ Assert (std::abs ((int) dataset_key.size () - (int) independent_values[0].size ()) < 2, ExcDataLostSync ());
+ }
// Look up the field name and get an
// iterator for the map. Doing this
// up front means that it only needs
// to check vector_name is in the map.
typename std::map <std::string, std::vector <std::vector <double> > >::iterator data_store_field = data_store.find(vector_name);
Assert (data_store_field != data_store.end(), ExcMessage("vector_name not in class"));
- // Repeat for field_components
- typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
- Assert (components != field_components.end(), ExcMessage("vector_name not in class"));
-
- Assert (components->second.size () <= dof_handler->get_fe ().n_components (), ExcIndexRange (components->second.size (), 0, dof_handler->get_fe ().n_components () + 1)); // 0 components legal as request to generate list internally!
- for (unsigned int component = 0; component < components->second.size (); component++)
- {
- Assert (components->second [component] < dof_handler->get_fe ().n_components (), ExcIndexRange (components->second[component], 0, dof_handler->get_fe ().n_components ()));
- } // hopefully the compiler will optimise
+ // Repeat for component_mask
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(vector_name);
+ Assert (mask != component_mask.end(), ExcMessage("vector_name not in class"));
+ Assert (mask->second.size () == dof_handler->get_fe ().n_components (), ExcDimensionMismatch (mask->second.size (), dof_handler->get_fe ().n_components ()));
+ unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true);
typename std::vector <internal::PointValueHistory::PointGeometryData <dim> >::iterator point = point_geometry_data.begin ();
Vector <double> value (dof_handler->get_fe().n_components());
// Make a Vector <double> for the value
// at the point. It will have as many
// components as there are in the fe.
- VectorTools::point_value (*dof_handler, solution, point->requested_location, value);
+ VectorTools::point_value (*dof_handler, solution, point->requested_location, value);
-
- // Look up the component to add
- // in field_components, and
- // access the data associated with
- // that component
- for (unsigned int comp = 0; comp < components->second.size(); comp++)
+ // Look up the component_mask and add
+ // in components according to that mask
+ for (unsigned int store_index = 0, comp = 0; comp < mask->second.size(); comp++)
{
- data_store_field->second[data_store_index * components->second.size() + comp].push_back (value (components->second[comp]));
+ if (mask->second[comp])
+ {
+ data_store_field->second[data_store_index * n_stored + store_index].push_back (value (comp));
+ store_index++;
+ }
}
}
}
for (; data_store_begin != data_store.end (); data_store_begin++)
{
- unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(data_store_begin->first);
+ unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true);
std::vector <std::string> names = (component_names_map.find (data_store_begin->first))->second;
+
if (names.size() > 0)
{
- for (unsigned int component = 0; component < names.size(); component++)
- {
- to_gnuplot << "<" << names[component] << "> ";
- }
+ AssertThrow (names.size() == n_stored, ExcDimensionMismatch (names.size(), n_stored));
+ for (unsigned int component = 0; component < names.size(); component++)
+ {
+ to_gnuplot << "<" << names[component] << "> ";
+ }
}
else
- {
- for (unsigned int component = 0; component < n_comps; component++)
+ {
+ for (unsigned int component = 0; component < n_stored; component++)
{
- to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
+ to_gnuplot << "<" << data_store_begin->first << "_" << component << "> ";
}
}
}
for (; data_store_begin != data_store.end (); data_store_begin++)
{
- unsigned int n_comps = (field_components.find (data_store_begin->first))->second.size(); // named to distinguish it from FE::n_components
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(data_store_begin->first);
+ unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true);
- for (unsigned int component = 0; component < n_comps; component++)
+ for (unsigned int component = 0; component < n_stored; component++)
{
- to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_comps + component][key];
+ to_gnuplot << " " << (data_store_begin->second)[data_store_index * n_stored + component][key];
}
}
to_gnuplot << "\n";
std::map <std::string, std::vector <std::vector <double> > >::iterator
data_store_begin = data_store.begin ();
+ if (data_store_begin != data_store.end())
+ {
+ out << "Mnemonic: data set size (mask size, n true components) : n data sets\n";
+ }
for (; data_store_begin != data_store.end (); data_store_begin++)
{
// Find field mnemonic
std::string vector_name = data_store_begin->first;
- typename std::map <std::string, std::vector <unsigned int> >::iterator components = field_components.find(vector_name);
- Assert (components != field_components.end(), ExcMessage("vector_name not in class"));
+ typename std::map <std::string, std::vector <bool> >::iterator mask = component_mask.find(vector_name);
+ Assert (mask != component_mask.end(), ExcMessage("vector_name not in class"));
typename std::map <std::string, std::vector <std::string> >::iterator component_names = component_names_map.find(vector_name);
Assert (component_names != component_names_map.end(), ExcMessage("vector_name not in class"));
if (data_store_begin->second.size () != 0)
{
out << data_store_begin->first << ": " << data_store_begin->second.size () << " (";
- out << components->second.size() << ") : ";
+ out << mask->second.size() << ", " << std::count (mask->second.begin(), mask->second.end(), true) << ") : ";
out << (data_store_begin->second)[0].size () << "\n";
}
else
{
out << data_store_begin->first << ": " << data_store_begin->second.size () << " (";
- out << components->second.size() << ") : ";
+ out << mask->second.size() << ", " << std::count (mask->second.begin(), mask->second.end(), true) << ") : ";
out << "No points added" << "\n";
}
// add names, if available
void PointValueHistory<dim>
::tria_change_listener ()
{
- // this function is called by the
- // Triangulation whenever something
- // changes, by virtue of having
- // attached the function to the
- // signal handler in the
- // triangulation object
-
- // we record the fact that the mesh
- // has changed. we need to take
- // this into account next time we
- // evaluate the solution
+ // this function is called by the
+ // Triangulation whenever something
+ // changes, by virtue of having
+ // attached the function to the
+ // signal handler in the
+ // triangulation object
+
+ // we record the fact that the mesh
+ // has changed. we need to take
+ // this into account next time we
+ // evaluate the solution
triangulation_changed = true;
}
triangulation.refine_global(2); // refine 2 times to make 5 nodes per side
// make a DOF handler, a model solution filled with ones and a flow vector
- FESystem<dim> finite_element(FE_Q<dim > (1 + 1), dim, FE_Q<dim > (1), 1);
- DoFHandler <dim > dof_handler(triangulation);
+ // FESystem<dim> finite_element(FE_Q<dim > (1 + 1), dim, FE_Q<dim > (1), 1); // Already done!
+ // DoFHandler <dim > dof_handler(triangulation);
dof_handler.distribute_dofs(finite_element);
DoFRenumbering::Cuthill_McKee(dof_handler);
std::vector < std::vector <Point <dim> > > selected_locations;
node_monitor.get_points(selected_locations);
- // write output to a file
+ // write output to a file
Vector<double> node_locations = node_monitor.mark_locations();
- // write output to a file
+ // write output to a file
}
double delta_t = 0.000001;
deallog << "Starting data files" << std::endl;
- // copy all the data into deallog and
- // delete those files
+ // copy all the data into deallog and
+ // delete those files
const std::string filenames[]
= { "point_value_history_01/node_00.gpl",
- "point_value_history_01/node_01.gpl",
- "point_value_history_01/node_02.gpl",
- "point_value_history_01/node_03.gpl",
- "point_value_history_01/node_04.gpl",
- "point_value_history_01/node_05.gpl",
- "point_value_history_01/node_indep.gpl",
- "point_value_history_01/Test_Copy_00.gpl",
- "point_value_history_01/Test_Copy_indep.gpl",
- "point_value_history_01/no_dof_indep.gpl" };
+ "point_value_history_01/node_01.gpl",
+ "point_value_history_01/node_02.gpl",
+ "point_value_history_01/node_03.gpl",
+ "point_value_history_01/node_04.gpl",
+ "point_value_history_01/node_05.gpl",
+ "point_value_history_01/node_indep.gpl",
+ "point_value_history_01/Test_Copy_00.gpl",
+ "point_value_history_01/Test_Copy_indep.gpl",
+ "point_value_history_01/no_dof_indep.gpl" };
for (unsigned int i=0; i<sizeof(filenames)/sizeof(filenames[0]); ++i)
- {
- deallog << "Copying output file " << filenames[i]
- << std::endl;
+ {
+ deallog << "Copying output file " << filenames[i]
+ << std::endl;
- std::ifstream in(filenames[i].c_str());
- AssertThrow (in, ExcIO());
+ std::ifstream in(filenames[i].c_str());
+ AssertThrow (in, ExcIO());
- std::string s;
- while (in)
- {
- std::getline (in, s);
- deallog << s << std::endl;
- }
+ std::string s;
+ while (in)
+ {
+ std::getline (in, s);
+ deallog << s << std::endl;
+ }
- std::remove (filenames[i].c_str());
+ std::remove (filenames[i].c_str());
- deallog << std::endl;
- }
+ deallog << std::endl;
+ }
}
Closed: 1
Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 1
Geometric Data
# Requested location: 1.0 0.20
# DoF_index : Support location (for each component)
Independent value(s): 1 : 1
-Solution: 3 : 1
+Mnemonic: data set size (mask size, n true components) : n data sets
+Solution: 3 (3, 3) : 1
***end of status output***
Closed: 1
Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 1
Geometric Data
# Requested location: 0.0 0.0
# DoF_index : Support location (for each component)
Independent value(s): 1 : 10
-Post Processed Vector: 18 : 10
-Solution: 18 : 10
+Mnemonic: data set size (mask size, n true components) : n data sets
+Post Processed Vector: 18 (3, 3) : 10
+Solution: 18 (3, 3) : 10
***end of status output***
Closed: 1
Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 0
Geometric Data
No points stored currently
--- /dev/null
+//---------------------------- point_value_history_02.cc ---------------------------
+// $Id: point_value_history_02.cc 23710 2011-05-17 04:50:10Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2009, 2010 by the deal.II authors and Michael Rapson
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- point_value_history_02.cc ---------------------------
+
+
+
+// ************************************************
+// A test program for the PointValueHistory class
+// Currently this only tests a finite element system
+// with 3 components, on a hyper cube and with Vectors.
+// The testing is done in two dimensions.
+// ************************************************
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+
+#include <deal.II/numerics/point_value_history.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+
+template <int dim>
+class Postprocess : public DataPostprocessor<dim>
+{
+ public:
+
+ void compute_derived_quantities_vector (
+ const std::vector< Vector< double > > &,
+ const std::vector< std::vector< Tensor< 1, dim > > > &,
+ const std::vector< std::vector< Tensor< 2, dim > > > &,
+ const std::vector< Point< dim > > &,
+ const std::vector< Point< dim > > &,
+ std::vector< Vector< double > > &
+ ) const;
+
+ std::vector<std::string> get_names () const;
+ UpdateFlags get_needed_update_flags () const;
+ unsigned int n_output_variables () const;
+ // The following function is not required
+ // by the point_value_history class.
+ //std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ // get_data_component_interpretation () const;
+};
+
+template <int dim>
+std::vector<std::string>
+Postprocess<dim>::get_names() const
+{
+ std::vector<std::string> names;
+ names.push_back ("Vector_out");
+ names.push_back ("Vector_out_1");
+ names.push_back ("Vector_out_2");
+ names.push_back ("Scalar_out");
+ return names;
+}
+
+template <int dim>
+UpdateFlags
+Postprocess<dim>::get_needed_update_flags () const
+{
+ return update_values | update_gradients | update_hessians;
+}
+
+template <int dim>
+unsigned int
+Postprocess<dim>::n_output_variables () const
+{
+ return 4;
+}
+
+
+template <int dim>
+void
+Postprocess<dim>::compute_derived_quantities_vector (
+ const std::vector< Vector< double > > & uh,
+ const std::vector< std::vector< Tensor< 1, dim > > > & duh,
+ const std::vector< std::vector< Tensor< 2, dim > > > & dduh,
+ const std::vector< Point< dim > > & /* normals */,
+ const std::vector< Point< dim > > & /* locations */,
+ std::vector< Vector< double > > & computed_quantities
+ ) const
+{
+ Assert(computed_quantities.size() == uh.size(),
+ ExcDimensionMismatch (computed_quantities.size(), uh.size()));
+
+ for (unsigned int i=0; i<computed_quantities.size(); i++)
+ {
+ Assert(computed_quantities[i].size() == 4,
+ ExcDimensionMismatch (computed_quantities[i].size(), 3));
+ Assert(uh[i].size() == dim + 1, ExcDimensionMismatch (uh[i].size(), dim + 1));
+
+ computed_quantities[i](0) = duh[i][0].norm(); // norm of x gradient
+ computed_quantities[i](1) = dduh[i][1].norm(); // norm of y hessian
+ computed_quantities[i](2) = duh[i][0].norm() + dduh[i][1].norm(); // norm of y hessian
+ computed_quantities[i](3) = uh[i].l2_norm();
+ }
+}
+
+
+
+template <int dim>
+class TestPointValueHistory
+{
+public:
+ TestPointValueHistory();
+ void run();
+
+private:
+ void output_results (unsigned int step, Vector <double> solution) const;
+
+ Triangulation <dim> triangulation;
+ FESystem<dim> finite_element;
+ DoFHandler <dim > dof_handler;
+ PointValueHistory <dim> test_copy;
+ std::vector <Point <dim> > postprocessor_locations;
+};
+
+
+
+
+template <int dim>
+TestPointValueHistory<dim>::TestPointValueHistory() :
+finite_element(FE_Q<dim > (1 + 1), dim, FE_Q<dim > (1), 1),
+
+dof_handler(triangulation)
+{ }
+
+
+
+
+template <int dim>
+void TestPointValueHistory<dim>::run()
+{
+ // Make a triangulation
+ GridGenerator::hyper_cube(triangulation, 0, 1);
+ triangulation.refine_global(2); // refine 2 times to make 5 nodes per side
+
+ // make a DOF handler, a model solution filled with ones and a flow vector
+// FESystem<dim> finite_element(FE_Q<dim > (1 + 1), dim, FE_Q<dim > (1), 1);
+// DoFHandler <dim > dof_handler(triangulation);
+ dof_handler.distribute_dofs(finite_element);
+ DoFRenumbering::Cuthill_McKee(dof_handler);
+
+ // renumber for components so that same dof indices are used for BlockVectors and normal Vectors
+ std::vector<unsigned int> block_component(dim + 1, 0);
+ block_component[dim] = 1; // component dim = pressure component!
+ DoFRenumbering::component_wise(dof_handler, block_component);
+
+ // Vector
+ Vector <double> solution, post_processed, poles;
+ solution.reinit(dof_handler.n_dofs());
+ post_processed.reinit(dof_handler.n_dofs());
+ poles.reinit(dof_handler.n_dofs());
+
+// // BlockVector
+// std::vector<unsigned int> dofs_per_block(2);
+// DoFTools::count_dofs_per_block(dof_handler, dofs_per_block, block_component);
+// const unsigned int n_u = dofs_per_block[0],
+// n_p = dofs_per_block[1];
+// BlockVector<double> solution;
+// solution.reinit(2);
+// solution.block(0).reinit(n_u);
+// solution.block(1).reinit(n_p);
+// solution.collect_sizes();
+//
+// BlockVector<double> post_processed;
+// post_processed.reinit(2);
+// post_processed.block(0).reinit(n_u);
+// post_processed.block(1).reinit(n_p);
+// post_processed.collect_sizes();
+// BlockVector<double> poles;
+// poles.reinit(2);
+// poles.block(0).reinit(n_u);
+// poles.block(1).reinit(n_p);
+// poles.collect_sizes();
+
+ // set up a simple linear discrete time system so that time plots vary
+ // over the mesh but can be easily checked. The basic idea is to have each
+ // component of the fe_system to depend on a specific dimension (i.e component 0
+ // depends on dim 0, etc. % dim handles the case where there are more components
+ // than dimensions. The code breaks down at the edges of the mesh and this is
+ // not corrected for. The code used in this test simplified from point_value_history_01.
+ {
+ Quadrature<dim> quadrature_formula(finite_element.get_unit_support_points ());
+ FEValues<dim> fe_values(finite_element, quadrature_formula, update_values | update_quadrature_points); // just need local_dof_indices and quadrature_points
+
+ std::vector<unsigned int> local_dof_indices(finite_element.dofs_per_cell);
+ std::vector<Point<dim> > dof_locations(finite_element.dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell, endc;
+ cell = dof_handler.begin_active();
+ endc = dof_handler.end();
+ for (; cell != endc; ++cell)
+ {
+ fe_values.reinit(cell); // need to get local_dof_indices
+ cell->get_dof_indices(local_dof_indices);
+ dof_locations = fe_values.get_quadrature_points();
+
+ for (unsigned int dof = 0; dof != finite_element.dofs_per_cell; dof++)
+ {
+ unsigned int dof_component = finite_element.system_to_component_index(dof).first;
+
+ poles(local_dof_indices[dof]) = -dof_locations[dof](dof_component % dim);
+
+ if (dof_component == dim) // components start numbering at 0
+ poles(local_dof_indices[dof]) = -0.1; // dim+1th component is not handled well by the code above
+
+ solution(local_dof_indices [dof]) = 1; // start all solutions at 1
+ }
+
+ } // loop over all cells
+ poles.add(1.0); // slow down the pole settling time
+ post_processed.add(3.0); // set to starting value.
+// output_results(10, poles);
+ }
+
+ // Setup monitor node to print variation over time
+ unsigned int n_inputs = 1;
+ PointValueHistory<dim> node_monitor (dof_handler, n_inputs);
+ PointValueHistory<dim> no_dof_handler (n_inputs);
+
+ // check that the assignment operator is valid
+ test_copy = node_monitor;
+ test_copy.add_point(Point<2>(1, 0.2));
+ test_copy.add_field_name("Solution");
+ std::vector < std::vector <Point <dim> > > selected_locations;
+ test_copy.get_support_locations(selected_locations);
+ test_copy.mark_support_locations();
+ test_copy.close();
+ test_copy.start_new_dataset(0.1);
+// triangulation.refine_global(1); // should mark the triangulation as changed
+ test_copy.evaluate_field("Solution", solution);
+ std::vector <double> input_value(n_inputs, 1);
+ test_copy.push_back_independent(input_value);
+ test_copy.write_gnuplot("point_value_history_02/Test_Copy");
+ test_copy.status(deallog.get_file_stream());
+ test_copy.clear ();
+ // end of assignment operator check
+
+ {
+ node_monitor.add_field_name("Solution");
+ std::vector <std::string> solution_names;
+ solution_names.push_back("X velocity");
+ solution_names.push_back("Y velocity");
+ solution_names.push_back("Z velocity");
+ node_monitor.add_component_names ("Solution", solution_names);
+ node_monitor.add_field_name("Post Processed Vector"); // not sensitive to spaces
+ std::vector <bool> component_mask (3, false);
+ component_mask[2] = true;
+ node_monitor.add_field_name("Pressure", component_mask);
+ component_mask = std::vector <bool> (3, false);
+ component_mask[1] = true;
+ node_monitor.add_field_name("Req_sol", component_mask);
+ component_mask = std::vector <bool> (4, true);
+ component_mask[3] = false;
+ node_monitor.add_field_name("Vector_out", component_mask);
+ component_mask = std::vector <bool> (4, false);
+ component_mask[3] = true;
+ node_monitor.add_field_name("Scalar_out", component_mask);
+
+ std::vector <std::string> indep_names;
+ indep_names.push_back ("Input");
+ node_monitor.add_independent_names(indep_names);
+
+ // two alternatives here, adding a point at a time or a vector of points
+ // 2d points
+ std::vector <Point < 2 > > point_vector(5, Point < 2 > ());
+ point_vector[0] = Point < 2 > (0, 0); // some of these points will hit a node, others won't
+ point_vector[1] = Point < 2 > (0.25, 0);
+ point_vector[2] = Point < 2 > (0.25, 0.45);
+ point_vector[3] = Point < 2 > (0.45, 0.45);
+ point_vector[4] = Point < 2 > (0.8, 0.8);
+
+ node_monitor.add_points(point_vector);
+ node_monitor.add_point(Point<2>(1, 0.2)); // add a single point
+
+ // MonitorNode requires that the instance is 'closed' before any data is added
+ // this ensures that points are not added once time starts.
+ node_monitor.close();
+ no_dof_handler.close(); // closing still required!
+
+ std::vector < std::vector <Point <dim> > > selected_locations;
+ node_monitor.get_support_locations(selected_locations);
+ Vector<double> node_locations = node_monitor.mark_support_locations();
+ QGauss<dim> postprocess_quadrature (2);
+ node_monitor.get_postprocessor_locations(postprocess_quadrature, postprocessor_locations);
+ }
+
+ double delta_t = 0.000001;
+ double t_max = 0.00001;
+ unsigned int step = 0;
+
+ for (double time = 0; time < t_max; time = time + delta_t)
+ {
+ node_monitor.start_new_dataset(time);
+ no_dof_handler.start_new_dataset(time);
+ // time and input are special, they don't vary over the mesh.
+
+ std::vector <double> input_value(n_inputs, 1); // manufacture an input value
+ node_monitor.push_back_independent(input_value);
+ no_dof_handler.push_back_independent(input_value);
+ node_monitor.evaluate_field("Solution", solution);
+ node_monitor.evaluate_field("Post Processed Vector", post_processed);
+ node_monitor.evaluate_field("Pressure", solution);
+ node_monitor.evaluate_field_at_requested_location("Req_sol", solution);
+
+ Postprocess<dim> postprocessor;
+ QGauss<dim> postprocess_quadrature (2);
+ std::vector<std::string> names;
+ names.push_back ("Vector_out");
+ names.push_back ("Scalar_out");
+ node_monitor.evaluate_field(names, solution, postprocessor, postprocess_quadrature);
+// output_results (step, solution);
+ step++;
+
+ solution.scale(poles); // decaying exponentials of varying time constants
+ post_processed = solution;
+ post_processed.add(2.0); // simple post processing, giving it a dc offset
+ }
+ triangulation.refine_global(1); // should mark the triangulation as changed
+ node_monitor.write_gnuplot("point_value_history_02/node", postprocessor_locations);
+ no_dof_handler.write_gnuplot("point_value_history_02/no_dof");
+
+ node_monitor.status (deallog.get_file_stream());
+ no_dof_handler.status (deallog.get_file_stream());
+
+ deallog << "Starting data files" << std::endl;
+
+ // copy all the data into deallog and
+ // delete those files
+ const std::string filenames[]
+ = { "point_value_history_02/node_00.gpl",
+ "point_value_history_02/node_01.gpl",
+ "point_value_history_02/node_02.gpl",
+ "point_value_history_02/node_03.gpl",
+ "point_value_history_02/node_04.gpl",
+ "point_value_history_02/node_05.gpl",
+ "point_value_history_02/node_indep.gpl",
+ "point_value_history_02/Test_Copy_00.gpl",
+ "point_value_history_02/Test_Copy_indep.gpl",
+ "point_value_history_02/no_dof_indep.gpl" };
+
+ for (unsigned int i=0; i<sizeof(filenames)/sizeof(filenames[0]); ++i)
+ {
+ deallog << "Copying output file " << filenames[i]
+ << std::endl;
+
+ std::ifstream in(filenames[i].c_str());
+ AssertThrow (in, ExcIO());
+
+ std::string s;
+ while (in)
+ {
+ std::getline (in, s);
+ deallog << s << std::endl;
+ }
+
+ std::remove (filenames[i].c_str());
+
+ deallog << std::endl;
+ }
+}
+
+template <int dim>
+void TestPointValueHistory<dim>::output_results (unsigned int step, Vector <double> solution) const
+{
+ std::vector<std::string> solution_names (dim, "velocity");
+ solution_names.push_back ("pressure");
+
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ data_component_interpretation
+ (dim, DataComponentInterpretation::component_is_part_of_vector);
+ data_component_interpretation
+ .push_back (DataComponentInterpretation::component_is_scalar);
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, solution_names,
+ DataOut<dim>::type_dof_data,
+ data_component_interpretation);
+ data_out.build_patches (2);
+
+ std::ostringstream filename;
+ filename << "point_value_history_02/solution-"
+ << Utilities::int_to_string (step, 2)
+ << ".gpl";
+
+ std::ofstream output (filename.str().c_str());
+ data_out.write_gnuplot (output);
+}
+
+
+
+int main()
+{
+ std::ofstream logfile("point_value_history_02/output");
+ logfile << std::setprecision(2);
+ deallog << std::setprecision(2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ TestPointValueHistory<2> test;
+ test.run();
+}
--- /dev/null
+
+***PointValueHistory status output***
+
+Closed: 1
+Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 1
+Geometric Data
+# Requested location: 1.0 0.20
+# DoF_index : Support location (for each component)
+138 : 1.0 0.25
+139 : 1.0 0.25
+181 : 1.0 0.25
+
+
+Independent value(s): 1 : 1
+Mnemonic: data set size (mask size, n true components) : n data sets
+Solution: 3 (3, 3) : 1
+
+***end of status output***
+
+***PointValueHistory status output***
+
+Closed: 1
+Cleared: 0
+Triangulation_changed: 1
+Have_dof_handler: 1
+Geometric Data
+# Requested location: 0.0 0.0
+# DoF_index : Support location (for each component)
+0 : 0.0 0.0
+1 : 0.0 0.0
+162 : 0.0 0.0
+
+# Requested location: 0.25 0.0
+# DoF_index : Support location (for each component)
+8 : 0.25 0.0
+9 : 0.25 0.0
+163 : 0.25 0.0
+
+# Requested location: 0.25 0.45
+# DoF_index : Support location (for each component)
+46 : 0.25 0.50
+47 : 0.25 0.50
+169 : 0.25 0.50
+
+# Requested location: 0.45 0.45
+# DoF_index : Support location (for each component)
+48 : 0.50 0.50
+49 : 0.50 0.50
+170 : 0.50 0.50
+
+# Requested location: 0.80 0.80
+# DoF_index : Support location (for each component)
+96 : 0.75 0.75
+97 : 0.75 0.75
+177 : 0.75 0.75
+
+# Requested location: 1.0 0.20
+# DoF_index : Support location (for each component)
+138 : 1.0 0.25
+139 : 1.0 0.25
+181 : 1.0 0.25
+
+
+Independent value(s): 1 : 10
+Names: <Input>
+Mnemonic: data set size (mask size, n true components) : n data sets
+Post Processed Vector: 18 (3, 3) : 10
+Pressure: 6 (3, 1) : 10
+Req_sol: 6 (3, 1) : 10
+Scalar_out: 6 (4, 1) : 10
+Solution: 18 (3, 3) : 10
+<X velocity> <Y velocity> <Z velocity>
+Vector_out: 18 (4, 3) : 10
+
+***end of status output***
+
+***PointValueHistory status output***
+
+Closed: 1
+Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 0
+Geometric Data
+No points stored currently
+
+Independent value(s): 1 : 10
+
+***end of status output***
+
+DEAL::Starting data files
+DEAL::Copying output file point_value_history_02/node_00.gpl
+DEAL::# Requested location: 0 0
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 0 : 0 0
+DEAL::# 1 : 0 0
+DEAL::# 162 : 0 0
+DEAL::# (Original components and locations, may be invalidated by mesh change.)
+DEAL::# Postprocessor location: 0.0528312 0.0528312 (may be approximate)
+DEAL::#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Post Processed Vector_1> <Post Processed Vector_2> <Pressure_0> <Req_sol_0> <Scalar_out_0> <X velocity> <Y velocity> <Z velocity> <Vector_out_0> <Vector_out_1> <Vector_out_2>
+DEAL::0 1 3 3 3 1 1 1.73205 1 1 1 0 1.19736e-08 1.19736e-08
+DEAL::1e-06 1 3 3 2.9 0.9 1 1.61377 1 1 0.9 1 1.16159e-08 1
+DEAL::2e-06 1 3 3 2.81 0.81 1 1.50525 1 1 0.81 1.89434 2 3.89434
+DEAL::3e-06 1 3 3 2.729 0.729 1 1.40644 1 1 0.729 2.69139 5.25 7.94139
+DEAL::4e-06 1 3 3 2.6561 0.6561 1 1.31707 1 1 0.6561 3.39968 9.21875 12.6184
+DEAL::5e-06 1 3 3 2.59049 0.59049 1 1.23666 1 1 0.59049 4.0276 13.5352 17.5628
+DEAL::6e-06 1 3 3 2.53144 0.531441 1 1.16468 1 1 0.531441 4.58314 17.9448 22.528
+DEAL::7e-06 1 3 3 2.4783 0.478297 1 1.10053 1 1 0.478297 5.07383 22.2779 27.3517
+DEAL::8e-06 1 3 3 2.43047 0.430467 1 1.04356 1 1 0.430467 5.50663 26.4253 31.9319
+DEAL::9e-06 1 3 3 2.38742 0.38742 1 0.99314 1 1 0.38742 5.88791 30.3212 36.2091
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/node_01.gpl
+DEAL::# Requested location: 0.25 0
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 8 : 0.25 0
+DEAL::# 9 : 0.25 0
+DEAL::# 163 : 0.25 0
+DEAL::# (Original components and locations, may be invalidated by mesh change.)
+DEAL::# Postprocessor location: 0.197169 0.0528312 (may be approximate)
+DEAL::#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Post Processed Vector_1> <Post Processed Vector_2> <Pressure_0> <Req_sol_0> <Scalar_out_0> <X velocity> <Y velocity> <Z velocity> <Vector_out_0> <Vector_out_1> <Vector_out_2>
+DEAL::0 1 3 3 3 1 1 1.73205 1 1 1 0 9.00988e-09 9.00988e-09
+DEAL::1e-06 1 2.75 3 2.9 0.9 1 1.53351 0.75 1 0.9 1 8.79108e-09 1
+DEAL::2e-06 1 2.5625 3 2.81 0.81 1 1.36981 0.5625 1 0.81 1.60566 2 3.60566
+DEAL::3e-06 1 2.42188 3 2.729 0.729 1 1.23359 0.421875 1 0.729 1.93361 5.25 7.18361
+DEAL::4e-06 1 2.31641 3 2.6561 0.6561 1 1.11938 0.316406 1 0.6561 2.06907 9.21875 11.2878
+DEAL::5e-06 1 2.2373 3 2.59049 0.59049 1 1.02303 0.237305 1 0.59049 2.07397 13.5352 15.6091
+DEAL::6e-06 1 2.17798 3 2.53144 0.531441 1 0.941352 0.177979 1 0.531441 1.99303 17.9448 19.9379
+DEAL::7e-06 1 2.13348 3 2.4783 0.478297 1 0.871889 0.133484 1 0.478297 1.8583 22.2779 24.1362
+DEAL::8e-06 1 2.10011 3 2.43047 0.430467 1 0.812685 0.100113 1 0.430467 1.69247 26.4253 28.1178
+DEAL::9e-06 1 2.07508 3 2.38742 0.38742 1 0.762164 0.0750847 1 0.38742 1.51142 30.3212 31.8326
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/node_02.gpl
+DEAL::# Requested location: 0.25 0.45
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 46 : 0.25 0.5
+DEAL::# 47 : 0.25 0.5
+DEAL::# 169 : 0.25 0.5
+DEAL::# (Original components and locations, may be invalidated by mesh change.)
+DEAL::# Postprocessor location: 0.197169 0.447169 (may be approximate)
+DEAL::#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Post Processed Vector_1> <Post Processed Vector_2> <Pressure_0> <Req_sol_0> <Scalar_out_0> <X velocity> <Y velocity> <Z velocity> <Vector_out_0> <Vector_out_1> <Vector_out_2>
+DEAL::0 1 3 3 3 1 1 1.73205 1 1 1 8.88178e-16 4.35771e-09 4.35771e-09
+DEAL::1e-06 1 2.75 2.5 2.9 0.9 0.55 1.32671 0.75 0.5 0.9 1 2.51837e-09 1
+DEAL::2e-06 1 2.5625 2.25 2.81 0.81 0.3025 1.07932 0.5625 0.25 0.81 1.60566 2 3.60566
+DEAL::3e-06 1 2.42188 2.125 2.729 0.729 0.165625 0.90924 0.421875 0.125 0.729 1.93361 3.75 5.68361
+DEAL::4e-06 1 2.31641 2.0625 2.6561 0.6561 0.0896875 0.780577 0.316406 0.0625 0.6561 2.06907 4.71875 6.78782
+DEAL::5e-06 1 2.2373 2.03125 2.59049 0.59049 0.0475586 0.677229 0.237305 0.03125 0.59049 2.07397 4.98047 7.05443
+DEAL::6e-06 1 2.17798 2.01562 2.53144 0.531441 0.0242896 0.591391 0.177979 0.015625 0.531441 1.99303 4.76123 6.75426
+DEAL::7e-06 1 2.13348 2.00781 2.4783 0.478297 0.0115738 0.51879 0.133484 0.0078125 0.478297 1.8583 4.2746 6.13289
+DEAL::8e-06 1 2.10011 2.00391 2.43047 0.430467 0.00476261 0.45676 0.100113 0.00390625 0.430467 1.69247 3.67699 5.36946
+DEAL::9e-06 1 2.07508 2.00195 2.38742 0.38742 0.00124056 0.40346 0.0750847 0.00195312 0.38742 1.51142 3.06777 4.57919
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/node_03.gpl
+DEAL::# Requested location: 0.45 0.45
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 48 : 0.5 0.5
+DEAL::# 49 : 0.5 0.5
+DEAL::# 170 : 0.5 0.5
+DEAL::# (Original components and locations, may be invalidated by mesh change.)
+DEAL::# Postprocessor location: 0.447169 0.447169 (may be approximate)
+DEAL::#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Post Processed Vector_1> <Post Processed Vector_2> <Pressure_0> <Req_sol_0> <Scalar_out_0> <X velocity> <Y velocity> <Z velocity> <Vector_out_0> <Vector_out_1> <Vector_out_2>
+DEAL::0 1 3 3 3 1 1 1.73205 1 1 1 8.88178e-16 4.35771e-09 4.35771e-09
+DEAL::1e-06 1 2.5 2.5 2.9 0.9 0.55 1.19216 0.5 0.5 0.9 1 2.51837e-09 1
+DEAL::2e-06 1 2.25 2.25 2.81 0.81 0.3025 0.918101 0.25 0.25 0.81 1.10566 2 3.10566
+DEAL::3e-06 1 2.125 2.125 2.729 0.729 0.165625 0.766829 0.125 0.125 0.729 0.916867 3.75 4.66687
+DEAL::4e-06 1 2.0625 2.0625 2.6561 0.6561 0.0896875 0.668761 0.0625 0.0625 0.6561 0.675079 4.71875 5.39383
+DEAL::5e-06 1 2.03125 2.03125 2.59049 0.59049 0.0475586 0.594518 0.03125 0.03125 0.59049 0.464784 4.98047 5.44525
+DEAL::6e-06 1 2.01562 2.01562 2.53144 0.531441 0.0242896 0.532629 0.015625 0.015625 0.531441 0.305802 4.76123 5.06703
+DEAL::7e-06 1 2.00781 2.00781 2.4783 0.478297 0.0115738 0.478603 0.0078125 0.0078125 0.478297 0.194193 4.2746 4.46879
+DEAL::8e-06 1 2.00391 2.00391 2.43047 0.430467 0.00476261 0.430527 0.00390625 0.00390625 0.430467 0.119462 3.67699 3.79646
+DEAL::9e-06 1 2.00195 2.00195 2.38742 0.38742 0.00124056 0.387426 0.00195312 0.00195312 0.38742 0.0711287 3.06777 3.1389
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/node_04.gpl
+DEAL::# Requested location: 0.8 0.8
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 96 : 0.75 0.75
+DEAL::# 97 : 0.75 0.75
+DEAL::# 177 : 0.75 0.75
+DEAL::# (Original components and locations, may be invalidated by mesh change.)
+DEAL::# Postprocessor location: 0.802831 0.802831 (may be approximate)
+DEAL::#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Post Processed Vector_1> <Post Processed Vector_2> <Pressure_0> <Req_sol_0> <Scalar_out_0> <X velocity> <Y velocity> <Z velocity> <Vector_out_0> <Vector_out_1> <Vector_out_2>
+DEAL::0 1 3 3 3 1 1 1.73205 1 1 1 0 1.19736e-08 1.19736e-08
+DEAL::1e-06 1 2.25 2.25 2.9 0.9 0.2 0.942205 0.25 0.25 0.9 1 2.64669e-09 1
+DEAL::2e-06 1 2.0625 2.0625 2.81 0.81 0.04 0.811864 0.0625 0.0625 0.81 0.394338 2 2.39434
+DEAL::3e-06 1 2.01562 2.01562 2.729 0.729 0.00875 0.729097 0.015625 0.015625 0.729 0.116627 0.75 0.866627
+DEAL::4e-06 1 2.00391 2.00391 2.6561 0.6561 0.00203125 0.656106 0.00390625 0.00390625 0.6561 0.0314119 0.21875 0.250162
+DEAL::5e-06 1 2.00098 2.00098 2.59049 0.59049 0.000488281 0.59049 0.000976562 0.000976562 0.59049 0.00813489 0.0585937 0.0667286
+DEAL::6e-06 1 2.00024 2.00024 2.53144 0.531441 0.000119629 0.531441 0.000244141 0.000244141 0.531441 0.00206896 0.0151367 0.0172057
+DEAL::7e-06 1 2.00006 2.00006 2.4783 0.478297 2.96021e-05 0.478297 6.10352e-05 6.10352e-05 0.478297 0.000521645 0.00384521 0.00436686
+DEAL::8e-06 1 2.00002 2.00002 2.43047 0.430467 7.36237e-06 0.430467 1.52588e-05 1.52588e-05 0.430467 0.000130962 0.000968933 0.00109989
+DEAL::9e-06 1 2 2 2.38742 0.38742 1.83582e-06 0.38742 3.8147e-06 3.8147e-06 0.38742 3.28093e-05 0.000243187 0.000275996
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/node_05.gpl
+DEAL::# Requested location: 1 0.2
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 138 : 1 0.25
+DEAL::# 139 : 1 0.25
+DEAL::# 181 : 1 0.25
+DEAL::# (Original components and locations, may be invalidated by mesh change.)
+DEAL::# Postprocessor location: 0.947169 0.197169 (may be approximate)
+DEAL::#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Post Processed Vector_1> <Post Processed Vector_2> <Pressure_0> <Req_sol_0> <Scalar_out_0> <X velocity> <Y velocity> <Z velocity> <Vector_out_0> <Vector_out_1> <Vector_out_2>
+DEAL::0 1 3 3 3 1 1 1.73205 1 1 1 8.88178e-16 4.35771e-09 4.35771e-09
+DEAL::1e-06 1 2 2.75 2.9 0.9 0.8 1.2072 0 0.75 0.9 1 3.606e-09 1
+DEAL::2e-06 1 2 2.5625 2.81 0.81 0.64 1.03515 0 0.5625 0.81 0.105662 2 2.10566
+DEAL::3e-06 1 2 2.42188 2.729 0.729 0.51125 0.893546 0 0.421875 0.729 0.00837341 5.25 5.25837
+DEAL::4e-06 1 2 2.31641 2.6561 0.6561 0.407031 0.775187 0 0.316406 0.6561 0.000161921 9.21875 9.21891
+DEAL::5e-06 1 2 2.2373 2.59049 0.59049 0.322168 0.675465 0 0.237305 0.59049 0.00032239 13.5352 13.5355
+DEAL::6e-06 1 2 2.17798 2.53144 0.531441 0.252659 0.590857 0 0.177979 0.531441 0.000115836 17.9448 17.9449
+DEAL::7e-06 1 2 2.13348 2.4783 0.478297 0.195398 0.518648 0 0.133484 0.478297 3.33639e-05 22.2779 22.2779
+DEAL::8e-06 1 2 2.10011 2.43047 0.430467 0.147964 0.456732 0 0.100113 0.430467 8.89157e-06 26.4253 26.4253
+DEAL::9e-06 1 2 2.07508 2.38742 0.38742 0.108462 0.403457 0 0.0750847 0.38742 2.29172e-06 30.3212 30.3212
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/node_indep.gpl
+DEAL::# Data independent of mesh location
+DEAL::# <Key> <Input>
+DEAL::0 1
+DEAL::1e-06 1
+DEAL::2e-06 1
+DEAL::3e-06 1
+DEAL::4e-06 1
+DEAL::5e-06 1
+DEAL::6e-06 1
+DEAL::7e-06 1
+DEAL::8e-06 1
+DEAL::9e-06 1
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/Test_Copy_00.gpl
+DEAL::# Requested location: 1 0.2
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 138 : 1 0.25
+DEAL::# 139 : 1 0.25
+DEAL::# 181 : 1 0.25
+DEAL::#
+DEAL::# <Key> <Indep_0> <Solution_0> <Solution_1> <Solution_2>
+DEAL::0.1 1 1 1 1
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/Test_Copy_indep.gpl
+DEAL::# Data independent of mesh location
+DEAL::# <Key> <Indep_0>
+DEAL::0.1 1
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_02/no_dof_indep.gpl
+DEAL::# Data independent of mesh location
+DEAL::# <Key> <Indep_0>
+DEAL::0 1
+DEAL::1e-06 1
+DEAL::2e-06 1
+DEAL::3e-06 1
+DEAL::4e-06 1
+DEAL::5e-06 1
+DEAL::6e-06 1
+DEAL::7e-06 1
+DEAL::8e-06 1
+DEAL::9e-06 1
+DEAL::
+DEAL::
--- /dev/null
+//---------------------------- point_value_history_03.cc ---------------------------
+// $Id: point_value_history_03.cc 23710 2011-05-17 04:50:10Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2009, 2010 by the deal.II authors and Michael Rapson
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- point_value_history_03.cc ---------------------------
+
+
+
+// ************************************************
+// A test program for the PointValueHistory class
+// Currently this only tests a finite element system
+// with 3 components, on a hyper cube and with Vectors.
+// The testing is done in two dimensions.
+// ************************************************
+
+
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/block_vector.h>
+
+#include <deal.II/numerics/point_value_history.h>
+
+#include <fstream>
+
+using namespace dealii;
+
+
+template <int dim>
+class Postprocess : public DataPostprocessor<dim>
+{
+ public:
+
+ void compute_derived_quantities_scalar (
+ const std::vector< double > &,
+ const std::vector< Tensor< 1, dim > > &,
+ const std::vector< Tensor< 2, dim > > &,
+ const std::vector< Point< dim > > &,
+ const std::vector< Point< dim > > &,
+ std::vector< Vector< double > > &
+ ) const;
+
+ std::vector<std::string> get_names () const;
+ UpdateFlags get_needed_update_flags () const;
+ unsigned int n_output_variables () const;
+ // The following function is not required
+ // by the point_value_history class.
+ //std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ // get_data_component_interpretation () const;
+};
+
+template <int dim>
+std::vector<std::string>
+Postprocess<dim>::get_names() const
+{
+ std::vector<std::string> names;
+ names.push_back ("X_gradient");
+ names.push_back ("X_hessian");
+ return names;
+}
+
+template <int dim>
+UpdateFlags
+Postprocess<dim>::get_needed_update_flags () const
+{
+ return update_values | update_gradients | update_hessians;
+}
+
+template <int dim>
+unsigned int
+Postprocess<dim>::n_output_variables () const
+{
+ return 2;
+}
+
+
+ template <int dim>
+ void
+ Postprocess<dim>::compute_derived_quantities_scalar (
+ const std::vector< double > & uh,
+ const std::vector< Tensor< 1, dim > > & duh,
+ const std::vector< Tensor< 2, dim > > & dduh,
+ const std::vector< Point< dim > > & /* normals */,
+ const std::vector< Point< dim > > & /* locations */,
+ std::vector< Vector< double > > & computed_quantities
+ ) const
+{
+ Assert(computed_quantities.size() == uh.size(),
+ ExcDimensionMismatch (computed_quantities.size(), uh.size()));
+
+ for (unsigned int i=0; i<computed_quantities.size(); i++)
+ {
+ Assert(computed_quantities[i].size() == 2,
+ ExcDimensionMismatch (computed_quantities[i].size(), 2));
+
+ computed_quantities[i](0) = duh[i][0]; // norm of x gradient
+ computed_quantities[i](1) = dduh[i][0].norm(); // norm of x hessian
+ }
+}
+
+
+
+template <int dim>
+class TestPointValueHistory
+{
+public:
+ TestPointValueHistory();
+ void run();
+
+private:
+ void output_results (unsigned int step, Vector <double> solution) const;
+
+ Triangulation <dim> triangulation;
+ FE_Q<dim> finite_element;
+ DoFHandler <dim > dof_handler;
+ PointValueHistory <dim> test_copy;
+ std::vector <Point <dim> > postprocessor_locations;
+};
+
+
+
+
+template <int dim>
+TestPointValueHistory<dim>::TestPointValueHistory() :
+ finite_element(2),
+ dof_handler(triangulation)
+{ }
+
+
+
+
+template <int dim>
+void TestPointValueHistory<dim>::run()
+{
+ // Make a triangulation
+ GridGenerator::hyper_cube(triangulation, 0, 1);
+ triangulation.refine_global(2); // refine 2 times to make 5 nodes per side
+
+ // make a DOF handler, a model solution filled with ones and a flow vector
+// DoFHandler <dim > dof_handler(triangulation);
+ dof_handler.distribute_dofs(finite_element);
+ DoFRenumbering::Cuthill_McKee(dof_handler);
+
+ // Vector
+ Vector <double> solution, post_processed, poles;
+ solution.reinit(dof_handler.n_dofs());
+ post_processed.reinit(dof_handler.n_dofs());
+ poles.reinit(dof_handler.n_dofs());
+
+// // BlockVector
+ // BlockVector not used with scalar fields generally
+
+ // set up a simple linear discrete time system so that time plots vary
+ // over the mesh but can be easily checked. The basic idea is to have each
+ // component of the fe_system to depend on a specific dimension (i.e component 0
+ // depends on dim 0, etc. % dim handles the case where there are more components
+ // than dimensions. The code breaks down at the edges of the mesh and this is
+ // not corrected for. The code used in this test is simplified from point_value_history_01.
+ {
+ Quadrature<dim> quadrature_formula(finite_element.get_unit_support_points ());
+ FEValues<dim> fe_values(finite_element, quadrature_formula, update_values | update_quadrature_points); // just need local_dof_indices and quadrature_points
+
+ std::vector<unsigned int> local_dof_indices(finite_element.dofs_per_cell);
+ std::vector<Point<dim> > dof_locations(finite_element.dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator cell, endc;
+ cell = dof_handler.begin_active();
+ endc = dof_handler.end();
+ for (; cell != endc; ++cell)
+ {
+ fe_values.reinit(cell); // need to get local_dof_indices
+ cell->get_dof_indices(local_dof_indices);
+ dof_locations = fe_values.get_quadrature_points();
+
+ for (unsigned int dof = 0; dof != finite_element.dofs_per_cell; dof++)
+ {
+ unsigned int dof_component = finite_element.system_to_component_index(dof).first;
+ // The line above will always evaluate to 1,
+ // simplifying the remaining lines
+
+ poles(local_dof_indices[dof]) = -dof_locations[dof](dof_component % dim);
+
+ if (dof_component == dim) // components start numbering at 0
+ poles(local_dof_indices[dof]) = -0.1; // dim+1th component is not handled well by the code above
+
+ solution(local_dof_indices [dof]) = 1; // start all solutions at 1
+ }
+ } // loop over all cells
+
+ poles.add(1.0); // slow down the pole settling time
+ post_processed.add(3.0); // set to starting value.
+// output_results(10, poles);
+ }
+
+ // Setup monitor node to print variation over time
+ unsigned int n_inputs = 1;
+ PointValueHistory<dim> node_monitor (dof_handler, n_inputs);
+ PointValueHistory<dim> no_dof_handler (n_inputs);
+
+ // check that the assignment operator is valid
+ test_copy = node_monitor;
+ test_copy.add_point(Point<2>(1, 0.2));
+ test_copy.add_field_name("Solution");
+ std::vector < std::vector <Point <dim> > > selected_locations;
+ test_copy.get_support_locations(selected_locations);
+ test_copy.mark_support_locations();
+ test_copy.close();
+ test_copy.start_new_dataset(0.1);
+ test_copy.evaluate_field("Solution", solution);
+ std::vector <double> input_value(n_inputs, 1);
+ test_copy.push_back_independent(input_value);
+ test_copy.write_gnuplot("point_value_history_03/Test_Copy");
+ test_copy.status(deallog.get_file_stream());
+ test_copy.clear ();
+ // end of assignment operator check
+
+ {
+ node_monitor.add_field_name("Solution");
+ std::vector <std::string> solution_names;
+ solution_names.push_back("Solution");
+ node_monitor.add_component_names ("Solution", solution_names);
+ node_monitor.add_field_name("Post Processed Vector"); // not sensitive to spaces
+ node_monitor.add_field_name("Pressure", 1);
+ std::vector <bool> component_mask (1, true);
+ node_monitor.add_field_name("Req_sol", component_mask);
+ component_mask = std::vector <bool> (2, false);
+ component_mask[0] = true;
+ node_monitor.add_field_name("X_gradient", component_mask);
+ component_mask = std::vector <bool> (2, false);
+ component_mask[1] = true;
+ node_monitor.add_field_name("X_hessian", component_mask);
+ std::vector <std::string> indep_names;
+ indep_names.push_back ("Input");
+ node_monitor.add_independent_names(indep_names);
+
+ // two alternatives here, adding a point at a time or a vector of points
+ // 2d points
+ std::vector <Point < 2 > > point_vector(5, Point < 2 > ());
+ point_vector[0] = Point < 2 > (0, 0); // some of these points will hit a node, others won't
+ point_vector[1] = Point < 2 > (0.25, 0);
+ point_vector[2] = Point < 2 > (0.25, 0.45);
+ point_vector[3] = Point < 2 > (0.45, 0.45);
+ point_vector[4] = Point < 2 > (0.8, 0.8);
+
+ node_monitor.add_points(point_vector);
+ node_monitor.add_point(Point<2>(1, 0.2)); // add a single point
+
+ // MonitorNode requires that the instance is 'closed' before any data is added
+ // this ensures that points are not added once time starts.
+ node_monitor.close();
+ no_dof_handler.close(); // closing still required!
+
+ std::vector < std::vector <Point <dim> > > selected_locations;
+ node_monitor.get_support_locations(selected_locations);
+ Vector<double> node_locations = node_monitor.mark_support_locations();
+ QGauss<dim> postprocess_quadrature (2);
+ node_monitor.get_postprocessor_locations(postprocess_quadrature, postprocessor_locations);
+ }
+
+ double delta_t = 0.000001;
+ double t_max = 0.00001;
+ unsigned int step = 0;
+
+ for (double time = 0; time < t_max; time = time + delta_t)
+ {
+ node_monitor.start_new_dataset(time);
+ no_dof_handler.start_new_dataset(time);
+ // time and input are special, they don't vary over the mesh.
+
+ std::vector <double> input_value(n_inputs, 1); // manufacture an input value
+ node_monitor.push_back_independent(input_value);
+ no_dof_handler.push_back_independent(input_value);
+ node_monitor.evaluate_field("Solution", solution);
+ node_monitor.evaluate_field("Post Processed Vector", post_processed);
+ node_monitor.evaluate_field("Pressure", solution);
+ node_monitor.evaluate_field_at_requested_location("Req_sol", solution);
+
+ Postprocess<dim> postprocessor;
+ QGauss<dim> postprocess_quadrature (2);
+ std::vector<std::string> names;
+ names.push_back ("X_gradient");
+ names.push_back ("X_hessian");
+ node_monitor.evaluate_field(names, solution, postprocessor, postprocess_quadrature);
+
+// output_results (step, solution);
+ step++;
+
+ solution.scale(poles); // decaying exponentials of varying time constants
+ post_processed = solution;
+ post_processed.add(2.0); // simple post processing, giving it a dc offset
+ }
+ node_monitor.write_gnuplot("point_value_history_03/node", postprocessor_locations);
+ no_dof_handler.write_gnuplot("point_value_history_03/no_dof"); // no point in adding postprocessor_locations
+
+ node_monitor.status (deallog.get_file_stream());
+ no_dof_handler.status (deallog.get_file_stream());
+
+ deallog << "Starting data files" << std::endl;
+
+ // copy all the data into deallog and
+ // delete those files
+ const std::string filenames[]
+ = { "point_value_history_03/node_00.gpl",
+ "point_value_history_03/node_01.gpl",
+ "point_value_history_03/node_02.gpl",
+ "point_value_history_03/node_03.gpl",
+ "point_value_history_03/node_04.gpl",
+ "point_value_history_03/node_05.gpl",
+ "point_value_history_03/node_indep.gpl",
+ "point_value_history_03/Test_Copy_00.gpl",
+ "point_value_history_03/Test_Copy_indep.gpl",
+ "point_value_history_03/no_dof_indep.gpl" };
+
+ for (unsigned int i=0; i<sizeof(filenames)/sizeof(filenames[0]); ++i)
+ {
+ deallog << "Copying output file " << filenames[i]
+ << std::endl;
+
+ std::ifstream in(filenames[i].c_str());
+ AssertThrow (in, ExcIO());
+
+ std::string s;
+ while (in)
+ {
+ std::getline (in, s);
+ deallog << s << std::endl;
+ }
+
+ std::remove (filenames[i].c_str());
+
+ deallog << std::endl;
+ }
+}
+
+template <int dim>
+void TestPointValueHistory<dim>::output_results (unsigned int step, Vector <double> solution) const
+{
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+
+ data_out.build_patches (2);
+
+ std::ostringstream filename;
+ filename << "point_value_history_03/solution-"
+ << Utilities::int_to_string (step, 2)
+ << ".gpl";
+
+ std::ofstream output (filename.str().c_str());
+ data_out.write_gnuplot (output);
+}
+
+
+
+int main()
+{
+ std::ofstream logfile("point_value_history_03/output");
+ logfile << std::setprecision(2);
+ deallog << std::setprecision(2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ TestPointValueHistory<2> test;
+ test.run();
+}
--- /dev/null
+
+***PointValueHistory status output***
+
+Closed: 1
+Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 1
+Geometric Data
+# Requested location: 1.0 0.20
+# DoF_index : Support location (for each component)
+69 : 1.0 0.25
+
+
+Independent value(s): 1 : 1
+Mnemonic: data set size (mask size, n true components) : n data sets
+Solution: 1 (1, 1) : 1
+
+***end of status output***
+
+***PointValueHistory status output***
+
+Closed: 1
+Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 1
+Geometric Data
+# Requested location: 0.0 0.0
+# DoF_index : Support location (for each component)
+0 : 0.0 0.0
+
+# Requested location: 0.25 0.0
+# DoF_index : Support location (for each component)
+4 : 0.25 0.0
+
+# Requested location: 0.25 0.45
+# DoF_index : Support location (for each component)
+23 : 0.25 0.50
+
+# Requested location: 0.45 0.45
+# DoF_index : Support location (for each component)
+24 : 0.50 0.50
+
+# Requested location: 0.80 0.80
+# DoF_index : Support location (for each component)
+48 : 0.75 0.75
+
+# Requested location: 1.0 0.20
+# DoF_index : Support location (for each component)
+69 : 1.0 0.25
+
+
+Independent value(s): 1 : 10
+Names: <Input>
+Mnemonic: data set size (mask size, n true components) : n data sets
+Post Processed Vector: 6 (1, 1) : 10
+Pressure: 6 (1, 1) : 10
+Req_sol: 6 (1, 1) : 10
+Solution: 6 (1, 1) : 10
+<Solution>
+X_gradient: 6 (2, 1) : 10
+X_hessian: 6 (2, 1) : 10
+
+***end of status output***
+
+***PointValueHistory status output***
+
+Closed: 1
+Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 0
+Geometric Data
+No points stored currently
+
+Independent value(s): 1 : 10
+
+***end of status output***
+
+DEAL::Starting data files
+DEAL::Copying output file point_value_history_03/node_00.gpl
+DEAL::# Requested location: 0 0
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 0 : 0 0
+DEAL::# Postprocessor location: 0.0528312 0.0528312#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Pressure_0> <Req_sol_0> <Solution> <X_gradient_0> <X_hessian_0>
+DEAL::0 1 3 1 1 1 0 8.46659e-09
+DEAL::1e-06 1 3 1 1 1 -1 8.35385e-09
+DEAL::2e-06 1 3 1 1 1 -1.89434 2
+DEAL::3e-06 1 3 1 1 1 -2.69139 5.25
+DEAL::4e-06 1 3 1 1 1 -3.39968 9.21875
+DEAL::5e-06 1 3 1 1 1 -4.0276 13.5352
+DEAL::6e-06 1 3 1 1 1 -4.58314 17.9448
+DEAL::7e-06 1 3 1 1 1 -5.07383 22.2779
+DEAL::8e-06 1 3 1 1 1 -5.50663 26.4253
+DEAL::9e-06 1 3 1 1 1 -5.88791 30.3212
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/node_01.gpl
+DEAL::# Requested location: 0.25 0
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 4 : 0.25 0
+DEAL::# Postprocessor location: 0.197169 0.0528312#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Pressure_0> <Req_sol_0> <Solution> <X_gradient_0> <X_hessian_0>
+DEAL::0 1 3 1 1 1 0 2.19707e-09
+DEAL::1e-06 1 2.75 0.75 0.75 0.75 -1 1.7506e-09
+DEAL::2e-06 1 2.5625 0.5625 0.5625 0.5625 -1.60566 2
+DEAL::3e-06 1 2.42188 0.421875 0.421875 0.421875 -1.93361 5.25
+DEAL::4e-06 1 2.31641 0.316406 0.316406 0.316406 -2.06907 9.21875
+DEAL::5e-06 1 2.2373 0.237305 0.237305 0.237305 -2.07397 13.5352
+DEAL::6e-06 1 2.17798 0.177979 0.177979 0.177979 -1.99303 17.9448
+DEAL::7e-06 1 2.13348 0.133484 0.133484 0.133484 -1.8583 22.2779
+DEAL::8e-06 1 2.10011 0.100113 0.100113 0.100113 -1.69247 26.4253
+DEAL::9e-06 1 2.07508 0.0750847 0.0750847 0.0750847 -1.51142 30.3212
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/node_02.gpl
+DEAL::# Requested location: 0.25 0.45
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 23 : 0.25 0.5
+DEAL::# Postprocessor location: 0.197169 0.447169#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Pressure_0> <Req_sol_0> <Solution> <X_gradient_0> <X_hessian_0>
+DEAL::0 1 3 1 1 1 -8.88178e-16 3.08136e-09
+DEAL::1e-06 1 2.75 0.75 0.75 0.75 -1 2.4926e-09
+DEAL::2e-06 1 2.5625 0.5625 0.5625 0.5625 -1.60566 2
+DEAL::3e-06 1 2.42188 0.421875 0.421875 0.421875 -1.93361 5.25
+DEAL::4e-06 1 2.31641 0.316406 0.316406 0.316406 -2.06907 9.21875
+DEAL::5e-06 1 2.2373 0.237305 0.237305 0.237305 -2.07397 13.5352
+DEAL::6e-06 1 2.17798 0.177979 0.177979 0.177979 -1.99303 17.9448
+DEAL::7e-06 1 2.13348 0.133484 0.133484 0.133484 -1.8583 22.2779
+DEAL::8e-06 1 2.10011 0.100113 0.100113 0.100113 -1.69247 26.4253
+DEAL::9e-06 1 2.07508 0.0750847 0.0750847 0.0750847 -1.51142 30.3212
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/node_03.gpl
+DEAL::# Requested location: 0.45 0.45
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 24 : 0.5 0.5
+DEAL::# Postprocessor location: 0.447169 0.447169#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Pressure_0> <Req_sol_0> <Solution> <X_gradient_0> <X_hessian_0>
+DEAL::0 1 3 1 1 1 -8.88178e-16 3.08136e-09
+DEAL::1e-06 1 2.5 0.5 0.55 0.5 -1 1.72237e-09
+DEAL::2e-06 1 2.25 0.25 0.3025 0.25 -1.10566 2
+DEAL::3e-06 1 2.125 0.125 0.165625 0.125 -0.916867 3.75
+DEAL::4e-06 1 2.0625 0.0625 0.0896875 0.0625 -0.675079 4.71875
+DEAL::5e-06 1 2.03125 0.03125 0.0475586 0.03125 -0.464784 4.98047
+DEAL::6e-06 1 2.01562 0.015625 0.0242896 0.015625 -0.305802 4.76123
+DEAL::7e-06 1 2.00781 0.0078125 0.0115738 0.0078125 -0.194193 4.2746
+DEAL::8e-06 1 2.00391 0.00390625 0.00476261 0.00390625 -0.119462 3.67699
+DEAL::9e-06 1 2.00195 0.00195312 0.00124056 0.00195312 -0.0711287 3.06777
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/node_04.gpl
+DEAL::# Requested location: 0.8 0.8
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 48 : 0.75 0.75
+DEAL::# Postprocessor location: 0.802831 0.802831#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Pressure_0> <Req_sol_0> <Solution> <X_gradient_0> <X_hessian_0>
+DEAL::0 1 3 1 1 1 0 8.46659e-09
+DEAL::1e-06 1 2.25 0.25 0.2 0.25 -1 2.0107e-09
+DEAL::2e-06 1 2.0625 0.0625 0.04 0.0625 -0.394338 2
+DEAL::3e-06 1 2.01562 0.015625 0.00875 0.015625 -0.116627 0.75
+DEAL::4e-06 1 2.00391 0.00390625 0.00203125 0.00390625 -0.0314119 0.21875
+DEAL::5e-06 1 2.00098 0.000976562 0.000488281 0.000976562 -0.00813489 0.0585937
+DEAL::6e-06 1 2.00024 0.000244141 0.000119629 0.000244141 -0.00206896 0.0151367
+DEAL::7e-06 1 2.00006 6.10352e-05 2.96021e-05 6.10352e-05 -0.000521645 0.00384521
+DEAL::8e-06 1 2.00002 1.52588e-05 7.36237e-06 1.52588e-05 -0.000130962 0.000968933
+DEAL::9e-06 1 2 3.8147e-06 1.83582e-06 3.8147e-06 -3.28093e-05 0.000243187
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/node_05.gpl
+DEAL::# Requested location: 1 0.2
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 69 : 1 0.25
+DEAL::# Postprocessor location: 0.947169 0.197169#
+DEAL::# <Key> <Input> <Post Processed Vector_0> <Pressure_0> <Req_sol_0> <Solution> <X_gradient_0> <X_hessian_0>
+DEAL::0 1 3 1 1 1 -8.88178e-16 3.08136e-09
+DEAL::1e-06 1 2 0 0 0 -1 1.84763e-10
+DEAL::2e-06 1 2 0 0 0 -0.105662 2
+DEAL::3e-06 1 2 0 0 0 -0.00837341 0.75
+DEAL::4e-06 1 2 0 0 0 0.000161921 0.21875
+DEAL::5e-06 1 2 0 0 0 0.00032239 0.0585938
+DEAL::6e-06 1 2 0 0 0 0.000115836 0.0151367
+DEAL::7e-06 1 2 0 0 0 3.33639e-05 0.00384521
+DEAL::8e-06 1 2 0 0 0 8.89157e-06 0.000968933
+DEAL::9e-06 1 2 0 0 0 2.29172e-06 0.000243187
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/node_indep.gpl
+DEAL::# Data independent of mesh location
+DEAL::# <Key> <Input>
+DEAL::0 1
+DEAL::1e-06 1
+DEAL::2e-06 1
+DEAL::3e-06 1
+DEAL::4e-06 1
+DEAL::5e-06 1
+DEAL::6e-06 1
+DEAL::7e-06 1
+DEAL::8e-06 1
+DEAL::9e-06 1
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/Test_Copy_00.gpl
+DEAL::# Requested location: 1 0.2
+DEAL::# DoF_index : Support location (for each component)
+DEAL::# 69 : 1 0.25
+DEAL::#
+DEAL::# <Key> <Indep_0> <Solution_0>
+DEAL::0.1 1 1
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/Test_Copy_indep.gpl
+DEAL::# Data independent of mesh location
+DEAL::# <Key> <Indep_0>
+DEAL::0.1 1
+DEAL::
+DEAL::
+DEAL::Copying output file point_value_history_03/no_dof_indep.gpl
+DEAL::# Data independent of mesh location
+DEAL::# <Key> <Indep_0>
+DEAL::0 1
+DEAL::1e-06 1
+DEAL::2e-06 1
+DEAL::3e-06 1
+DEAL::4e-06 1
+DEAL::5e-06 1
+DEAL::6e-06 1
+DEAL::7e-06 1
+DEAL::8e-06 1
+DEAL::9e-06 1
+DEAL::
+DEAL::
--- /dev/null
+//--------------------------------------------------------------------
+// $Id: cochain_01.cc 25686 2012-07-04 14:18:43Z bangerth $
+//
+// Copyright (C) 2012 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//--------------------------------------------------------------------
+
+// Test the functions in integrators/laplace.h
+// Output matrices and assert consistency of residuals
+
+#include "../tests.h"
+#include "../lib/test_grids.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/integrators/laplace.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nedelec.h>
+
+using namespace LocalIntegrators::Laplace;
+
+template <int dim>
+void test_cell(const FEValuesBase<dim>& fev)
+{
+ const unsigned int n = fev.dofs_per_cell;
+ unsigned int d=fev.get_fe().n_components();
+ FullMatrix<double> M(n,n);
+ cell_matrix(M,fev);
+ M.print(deallog,8);
+
+ Vector<double> u(n), v(n), w(n);
+ std::vector<std::vector<Tensor<1,dim> > >
+ ugrad(d,std::vector<Tensor<1,dim> >(fev.n_quadrature_points));
+
+ std::vector<unsigned int> indices(n);
+ for (unsigned int i=0;i<n;++i)
+ indices[i] = i;
+
+ deallog << "Residuals";
+ for (unsigned int i=0;i<n;++i)
+ {
+ u = 0.;
+ u(i) = 1.;
+ w = 0.;
+ fev.get_function_gradients(u, indices, ugrad, true);
+ cell_residual(w, fev, make_slice(ugrad));
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << ' ' << w.l2_norm();
+ if (d==1)
+ {
+ cell_residual(w, fev, ugrad[0]);
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << " e" << w.l2_norm();
+ }
+ }
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void test_boundary(const FEValuesBase<dim>& fev)
+{
+ const unsigned int n = fev.dofs_per_cell;
+ unsigned int d=fev.get_fe().n_components();
+ FullMatrix<double> M(n,n);
+ nitsche_matrix(M, fev, 17);
+ M.print(deallog,8);
+
+ Vector<double> u(n), v(n), w(n);
+ std::vector<std::vector<double> >
+ uval (d,std::vector<double>(fev.n_quadrature_points)),
+ null_val(d,std::vector<double>(fev.n_quadrature_points, 0.));
+ std::vector<std::vector<Tensor<1,dim> > >
+ ugrad (d,std::vector<Tensor<1,dim> >(fev.n_quadrature_points));
+
+ std::vector<unsigned int> indices(n);
+ for (unsigned int i=0;i<n;++i)
+ indices[i] = i;
+
+ deallog << "Residuals";
+ for (unsigned int i=0;i<n;++i)
+ {
+ u = 0.;
+ u(i) = 1.;
+ w = 0.;
+ fev.get_function_values(u, indices, uval, true);
+ fev.get_function_gradients(u, indices, ugrad, true);
+ nitsche_residual(w, fev, make_slice(uval), make_slice(ugrad), make_slice(null_val), 17);
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << ' ' << w.l2_norm();
+ if (d==1)
+ {
+ nitsche_residual(w, fev, uval[0], ugrad[0], null_val[0], 17);
+ M.vmult(v,u);
+ w.add(-1., v);
+ deallog << " e" << w.l2_norm();
+ }
+ }
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void test_face(const FEValuesBase<dim>& fev1,
+ const FEValuesBase<dim>& fev2)
+{
+ const unsigned int n1 = fev1.dofs_per_cell;
+ const unsigned int n2 = fev1.dofs_per_cell;
+ unsigned int d=fev1.get_fe().n_components();
+ FullMatrix<double> M11(n1,n1);
+ FullMatrix<double> M12(n1,n2);
+ FullMatrix<double> M21(n2,n1);
+ FullMatrix<double> M22(n2,n2);
+
+ ip_matrix(M11, M12, M21, M22, fev1, fev2, 17);
+ deallog << "M11" << std::endl;
+ M11.print(deallog,8);
+ deallog << "M22" << std::endl;
+ M22.print(deallog,8);
+ deallog << "M12" << std::endl;
+ M12.print(deallog,8);
+ deallog << "M21" << std::endl;
+ M21.print(deallog,8);
+
+ Vector<double> u1(n1), v1(n1), w1(n1);
+ Vector<double> u2(n2), v2(n2), w2(n2);
+ std::vector<std::vector<double> >
+ u1val (d,std::vector<double>(fev1.n_quadrature_points)),
+ u2val (d,std::vector<double>(fev2.n_quadrature_points));
+ std::vector<std::vector<Tensor<1,dim> > >
+ u1grad(d,std::vector<Tensor<1,dim> >(fev1.n_quadrature_points)),
+ u2grad(d,std::vector<Tensor<1,dim> >(fev2.n_quadrature_points));
+
+ std::vector<unsigned int> indices1(n1), indices2(n2);
+ for (unsigned int i=0;i<n1;++i) indices1[i] = i;
+ for (unsigned int i=0;i<n2;++i) indices2[i] = i;
+
+ deallog << "Residuals"; for (unsigned int i=0;i<n1;++i) indices1[i] = i;
+
+ for (unsigned int i1=0;i1<n1;++i1)
+ for (unsigned int i2=0;i2<n2;++i2)
+ {
+ u1 = 0.;
+ u1(i1) = 1.;
+ w1 = 0.;
+ fev1.get_function_values(u1, indices1, u1val, true);
+ fev1.get_function_gradients(u1, indices1, u1grad, true);
+ u2 = 0.;
+ u2(i2) = 1.;
+ w2 = 0.;
+ fev2.get_function_values(u2, indices2, u2val, true);
+ fev2.get_function_gradients(u2, indices2, u2grad, true);
+ ip_residual(w1, w2, fev1, fev2,
+ make_slice(u1val), make_slice(u1grad),
+ make_slice(u2val), make_slice(u2grad),
+ 17);
+ M11.vmult(v1,u1); w1.add(-1., v1);
+ M12.vmult(v1,u2); w1.add(-1., v1);
+ M21.vmult(v2,u1); w2.add(-1., v2);
+ M22.vmult(v2,u2); w2.add(-1., v2);
+ deallog << ' ' << w1.l2_norm() + w2.l2_norm();
+ if (d==1)
+ {
+ ip_residual(w1, w2, fev1, fev2, u1val[0], u1grad[0], u2val[0], u2grad[0], 17);
+ M11.vmult(v1,u1); w1.add(-1., v1);
+ M12.vmult(v1,u2); w1.add(-1., v1);
+ M21.vmult(v2,u1); w2.add(-1., v2);
+ M22.vmult(v2,u2); w2.add(-1., v2);
+ deallog << " e" << w1.l2_norm() + w2.l2_norm();
+ }
+ }
+ deallog << std::endl;
+}
+
+
+template <int dim>
+void
+test_fe(Triangulation<dim>& tr, FiniteElement<dim>& fe)
+{
+ deallog << fe.get_name() << std::endl << "cell matrix" << std::endl;
+ QGauss<dim> quadrature(fe.tensor_degree()+1);
+ FEValues<dim> fev(fe, quadrature, update_gradients);
+
+ typename Triangulation<dim>::cell_iterator cell1 = tr.begin(1);
+ fev.reinit(cell1);
+ test_cell(fev);
+
+ QGauss<dim-1> face_quadrature(fe.tensor_degree()+1);
+ FEFaceValues<dim> fef1(fe, face_quadrature, update_values | update_gradients | update_normal_vectors);
+ for (unsigned int i=0;i<GeometryInfo<dim>::faces_per_cell;++i)
+ {
+ deallog << "boundary_matrix " << i << std::endl;
+ fef1.reinit(cell1, i);
+ test_boundary(fef1);
+ }
+
+ FEFaceValues<dim> fef2(fe, face_quadrature, update_values | update_gradients);
+ typename Triangulation<dim>::cell_iterator cell2 = cell1->neighbor(1);
+
+ deallog << "face_matrix " << 0 << std::endl;
+ cell1 = tr.begin(1);
+ fef1.reinit(cell1, 1);
+ fef2.reinit(cell2, 0);
+ test_face(fef1, fef2);
+}
+
+
+template <int dim>
+void
+test(Triangulation<dim>& tr)
+{
+ FE_DGQ<dim> q1(1);
+ test_fe(tr, q1);
+
+// FE_DGQ<dim> q2(2);
+// test_fe(tr, q2);
+
+ FE_Nedelec<dim> n1(1);
+ test_fe(tr, n1);
+}
+
+
+int main()
+{
+ const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Triangulation<2> tr2;
+ TestGrids::hypercube(tr2, 1);
+ test(tr2);
+
+// Triangulation<3> tr3;
+// TestGrids::hypercube(tr3, 1);
+// test(tr3);
+
+}
--- /dev/null
+
+DEAL::Triangulation hypercube 2D refinement 1 steps 4 active cells 5 total cells
+DEAL::FE_DGQ<2>(1)
+DEAL::cell matrix
+DEAL::0.67 -0.17 -0.17 -0.33
+DEAL::-0.17 0.67 -0.33 -0.17
+DEAL::-0.17 -0.33 0.67 -0.17
+DEAL::-0.33 -0.17 -0.17 0.67
+DEAL::Residuals 0 e0 0 e0 0 e0 0 e0
+DEAL::boundary_matrix 0
+DEAL::11. 0.33 5.3 0.17
+DEAL::0.33 0 0.17 0
+DEAL::5.3 0.17 11. 0.33
+DEAL::0.17 0 0.33 0
+DEAL::Residuals 0 e0 0 e0 0 e0 0 e0
+DEAL::boundary_matrix 1
+DEAL::0 0.33 0 0.17
+DEAL::0.33 11. 0.17 5.3
+DEAL::0 0.17 0 0.33
+DEAL::0.17 5.3 0.33 11.
+DEAL::Residuals 0 e0 0 e0 0 e0 0 e0
+DEAL::boundary_matrix 2
+DEAL::11. 5.3 0.33 0.17
+DEAL::5.3 11. 0.17 0.33
+DEAL::0.33 0.17 0 0
+DEAL::0.17 0.33 0 0
+DEAL::Residuals 0 e0 0 e0 0 e0 0 e0
+DEAL::boundary_matrix 3
+DEAL::0 0 0.33 0.17
+DEAL::0 0 0.17 0.33
+DEAL::0.33 0.17 11. 5.3
+DEAL::0.17 0.33 5.3 11.
+DEAL::Residuals 0 e0 0 e0 0 e0 0 e0
+DEAL::face_matrix 0
+DEAL::M11
+DEAL::0 0.17 0 0.083
+DEAL::0.17 5.3 0.083 2.7
+DEAL::0 0.083 0 0.17
+DEAL::0.083 2.7 0.17 5.3
+DEAL::M22
+DEAL::5.3 0.17 2.7 0.083
+DEAL::0.17 0 0.083 0
+DEAL::2.7 0.083 5.3 0.17
+DEAL::0.083 0 0.17 0
+DEAL::M12
+DEAL::-0.17 0 -0.083 0
+DEAL::-5.3 -0.17 -2.7 -0.083
+DEAL::-0.083 0 -0.17 0
+DEAL::-2.7 -0.083 -5.3 -0.17
+DEAL::M21
+DEAL::-0.17 -5.3 -0.083 -2.7
+DEAL::0 -0.17 0 -0.083
+DEAL::-0.083 -2.7 -0.17 -5.3
+DEAL::0 -0.083 0 -0.17
+DEAL::Residuals 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0 0 e0
+DEAL::FE_Nedelec<2>(1)
+DEAL::cell matrix
+DEAL::1.0 0 -1.0 0 0 0 0 0 0 0 0 0
+DEAL::0 5.0 0 1.0 0 0 0 0 0 0 0 -1.7
+DEAL::-1.0 0 1.0 0 0 0 0 0 0 0 0 0
+DEAL::0 1.0 0 5.0 0 0 0 0 0 0 0 -1.7
+DEAL::0 0 0 0 1.0 0 -1.0 0 0 0 0 0
+DEAL::0 0 0 0 0 5.0 0 1.0 0 -1.7 0 0
+DEAL::0 0 0 0 -1.0 0 1.0 0 0 0 0 0
+DEAL::0 0 0 0 0 1.0 0 5.0 0 -1.7 0 0
+DEAL::0 0 0 0 0 0 0 0 1.0 0 0 0
+DEAL::0 0 0 0 0 -1.7 0 -1.7 0 2.2 0 0
+DEAL::0 0 0 0 0 0 0 0 0 0 1.0 0
+DEAL::0 -1.7 0 -1.7 0 0 0 0 0 0 0 2.2
+DEAL::Residuals 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::boundary_matrix 0
+DEAL::32. 0 1.0 0 0 0 0 0 0 0 -1.7 0
+DEAL::0 32. 0 1.0 0 0 0 0 0 0 0 -1.7
+DEAL::1.0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 1.0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 0 11. -18. 5.7 -9.2 -4.9 8.0 0 0
+DEAL::0 0 0 0 -18. 30. -9.2 15. 8.0 -13. 0 0
+DEAL::0 0 0 0 5.7 -9.2 11. -18. -4.9 8.0 0 0
+DEAL::0 0 0 0 -9.2 15. -18. 30. 8.0 -13. 0 0
+DEAL::0 0 0 0 -4.9 8.0 -4.9 8.0 3.4 -5.5 0 0
+DEAL::0 0 0 0 8.0 -13. 8.0 -13. -5.5 9.0 0 0
+DEAL::-1.7 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 -1.7 0 0 0 0 0 0 0 0 0 0
+DEAL::Residuals 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::boundary_matrix 1
+DEAL::0 0 1.0 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 1.0 0 0 0 0 0 0 0 0
+DEAL::1.0 0 32. 0 0 0 0 0 0 0 -1.7 0
+DEAL::0 1.0 0 32. 0 0 0 0 0 0 0 -1.7
+DEAL::0 0 0 0 11. 18. 5.7 9.2 -4.9 -8.0 0 0
+DEAL::0 0 0 0 18. 30. 9.2 15. -8.0 -13. 0 0
+DEAL::0 0 0 0 5.7 9.2 11. 18. -4.9 -8.0 0 0
+DEAL::0 0 0 0 9.2 15. 18. 30. -8.0 -13. 0 0
+DEAL::0 0 0 0 -4.9 -8.0 -4.9 -8.0 3.4 5.5 0 0
+DEAL::0 0 0 0 -8.0 -13. -8.0 -13. 5.5 9.0 0 0
+DEAL::0 0 -1.7 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 -1.7 0 0 0 0 0 0 0 0
+DEAL::Residuals 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::boundary_matrix 2
+DEAL::11. -18. 5.7 -9.2 0 0 0 0 0 0 -4.9 8.0
+DEAL::-18. 30. -9.2 15. 0 0 0 0 0 0 8.0 -13.
+DEAL::5.7 -9.2 11. -18. 0 0 0 0 0 0 -4.9 8.0
+DEAL::-9.2 15. -18. 30. 0 0 0 0 0 0 8.0 -13.
+DEAL::0 0 0 0 32. 0 1.0 0 -1.7 0 0 0
+DEAL::0 0 0 0 0 32. 0 1.0 0 -1.7 0 0
+DEAL::0 0 0 0 1.0 0 0 0 0 0 0 0
+DEAL::0 0 0 0 0 1.0 0 0 0 0 0 0
+DEAL::0 0 0 0 -1.7 0 0 0 0 0 0 0
+DEAL::0 0 0 0 0 -1.7 0 0 0 0 0 0
+DEAL::-4.9 8.0 -4.9 8.0 0 0 0 0 0 0 3.4 -5.5
+DEAL::8.0 -13. 8.0 -13. 0 0 0 0 0 0 -5.5 9.0
+DEAL::Residuals 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::boundary_matrix 3
+DEAL::11. 18. 5.7 9.2 0 0 0 0 0 0 -4.9 -8.0
+DEAL::18. 30. 9.2 15. 0 0 0 0 0 0 -8.0 -13.
+DEAL::5.7 9.2 11. 18. 0 0 0 0 0 0 -4.9 -8.0
+DEAL::9.2 15. 18. 30. 0 0 0 0 0 0 -8.0 -13.
+DEAL::0 0 0 0 0 0 1.0 0 0 0 0 0
+DEAL::0 0 0 0 0 0 0 1.0 0 0 0 0
+DEAL::0 0 0 0 1.0 0 32. 0 -1.7 0 0 0
+DEAL::0 0 0 0 0 1.0 0 32. 0 -1.7 0 0
+DEAL::0 0 0 0 0 0 -1.7 0 0 0 0 0
+DEAL::0 0 0 0 0 0 0 -1.7 0 0 0 0
+DEAL::-4.9 -8.0 -4.9 -8.0 0 0 0 0 0 0 3.4 5.5
+DEAL::-8.0 -13. -8.0 -13. 0 0 0 0 0 0 5.5 9.0
+DEAL::Residuals 0 0 0 0 0 0 0 0 0 0 0 0
+DEAL::face_matrix 0
+DEAL::M11
+DEAL::0 0 0.50 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 0.50 0 0 0 0 0 0 0 0
+DEAL::0.50 0 16. 0 0 0 0 0 0 0 -0.87 0
+DEAL::0 0.50 0 16. 0 0 0 0 0 0 0 -0.87
+DEAL::0 0 0 0 5.7 9.2 2.8 4.6 -2.5 -4.0 0 0
+DEAL::0 0 0 0 9.2 15. 4.6 7.5 -4.0 -6.5 0 0
+DEAL::0 0 0 0 2.8 4.6 5.7 9.2 -2.5 -4.0 0 0
+DEAL::0 0 0 0 4.6 7.5 9.2 15. -4.0 -6.5 0 0
+DEAL::0 0 0 0 -2.5 -4.0 -2.5 -4.0 1.7 2.8 0 0
+DEAL::0 0 0 0 -4.0 -6.5 -4.0 -6.5 2.8 4.5 0 0
+DEAL::0 0 -0.87 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 -0.87 0 0 0 0 0 0 0 0
+DEAL::M22
+DEAL::16. 0 0.50 0 0 0 0 0 0 0 -0.87 0
+DEAL::0 16. 0 0.50 0 0 0 0 0 0 0 -0.87
+DEAL::0.50 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 0.50 0 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 0 5.7 -9.2 2.8 -4.6 -2.5 4.0 0 0
+DEAL::0 0 0 0 -9.2 15. -4.6 7.5 4.0 -6.5 0 0
+DEAL::0 0 0 0 2.8 -4.6 5.7 -9.2 -2.5 4.0 0 0
+DEAL::0 0 0 0 -4.6 7.5 -9.2 15. 4.0 -6.5 0 0
+DEAL::0 0 0 0 -2.5 4.0 -2.5 4.0 1.7 -2.8 0 0
+DEAL::0 0 0 0 4.0 -6.5 4.0 -6.5 -2.8 4.5 0 0
+DEAL::-0.87 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 -0.87 0 0 0 0 0 0 0 0 0 0
+DEAL::M12
+DEAL::-0.50 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 -0.50 0 0 0 0 0 0 0 0 0 0
+DEAL::-16. 0 -0.50 0 0 0 0 0 0 0 0.87 0
+DEAL::0 -16. 0 -0.50 0 0 0 0 0 0 0 0.87
+DEAL::0 0 0 0 -5.7 9.2 -2.8 4.6 2.5 -4.0 0 0
+DEAL::0 0 0 0 -9.2 15. -4.6 7.5 4.0 -6.5 0 0
+DEAL::0 0 0 0 -2.8 4.6 -5.7 9.2 2.5 -4.0 0 0
+DEAL::0 0 0 0 -4.6 7.5 -9.2 15. 4.0 -6.5 0 0
+DEAL::0 0 0 0 2.5 -4.0 2.5 -4.0 -1.7 2.8 0 0
+DEAL::0 0 0 0 4.0 -6.5 4.0 -6.5 -2.8 4.5 0 0
+DEAL::0.87 0 0 0 0 0 0 0 0 0 0 0
+DEAL::0 0.87 0 0 0 0 0 0 0 0 0 0
+DEAL::M21
+DEAL::-0.50 0 -16. 0 0 0 0 0 0 0 0.87 0
+DEAL::0 -0.50 0 -16. 0 0 0 0 0 0 0 0.87
+DEAL::0 0 -0.50 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 -0.50 0 0 0 0 0 0 0 0
+DEAL::0 0 0 0 -5.7 -9.2 -2.8 -4.6 2.5 4.0 0 0
+DEAL::0 0 0 0 9.2 15. 4.6 7.5 -4.0 -6.5 0 0
+DEAL::0 0 0 0 -2.8 -4.6 -5.7 -9.2 2.5 4.0 0 0
+DEAL::0 0 0 0 4.6 7.5 9.2 15. -4.0 -6.5 0 0
+DEAL::0 0 0 0 2.5 4.0 2.5 4.0 -1.7 -2.8 0 0
+DEAL::0 0 0 0 -4.0 -6.5 -4.0 -6.5 2.8 4.5 0 0
+DEAL::0 0 0.87 0 0 0 0 0 0 0 0 0
+DEAL::0 0 0 0.87 0 0 0 0 0 0 0 0
+DEAL::Residuals 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009, 2010 by the deal.II authors
+// Copyright (C) 2009, 2010, 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
parallel::distributed::Triangulation<dim> tr(MPI_COMM_WORLD);
-
+
GridGenerator::hyper_cube(tr);
-
+
tr.refine_global (1);
for (unsigned int step=0; step<5;++step)
{
dofh.distribute_dofs (fe);
IndexSet owned_set = dofh.locally_owned_dofs();
-
+
IndexSet dof_set;
DoFTools::extract_locally_active_dofs (dofh, dof_set);
TrilinosWrappers::MPI::Vector x_rel;
x_rel.reinit(relevant_set, MPI_COMM_WORLD);
- x_rel = 0;
x_rel.compress();
-
+
ConstraintMatrix cm(relevant_set);
DoFTools::make_hanging_node_constraints (dofh, cm);
std::vector<bool> velocity_mask (dim+1, true);
-
+
velocity_mask[dim] = false;
VectorTools::interpolate_boundary_values (dofh,
ZeroFunction<dim>(dim+1),
cm,
velocity_mask);
-
+
cm.close ();
TrilinosWrappers::MPI::Vector x_test;
x_test.reinit(x_rel);
-
+
x_test=x;
bool throwing=false;
try
- {
+ {
cm.distribute(x_test);
}
catch (ExcMessage e)
//l2_norm() not possible for ghosted vectors...
//double a=0;//x_test.l2_norm();
//double b=0;//x_rel.l2_norm();
-
+
/* if (myid==0)
deallog << a << " vs " << b << std::endl;
*/
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009, 2010 by the deal.II authors
+// Copyright (C) 2009, 2010, 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
TrilinosWrappers::MPI::Vector x_rel;
x_rel.reinit(relevant_set, MPI_COMM_WORLD);
- x_rel = 0;
x_rel.compress();
ConstraintMatrix cm(relevant_set);
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009, 2010 by the deal.II authors
+// Copyright (C) 2009, 2010, 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
trans.interpolate(x);
x_rel.reinit(relevant_set, MPI_COMM_WORLD);
- x_rel = 0;
x_rel.compress();
ConstraintMatrix cm(relevant_set);
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009, 2010, 2011 by the deal.II authors
+// Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
tr.execute_coarsening_and_refinement ();
}
-
+
DoFHandler<dim> dofh(tr);
static const FE_Q<dim> fe(1);
dofh.distribute_dofs (fe);
IndexSet owned_set = dofh.locally_owned_dofs();
-
+
IndexSet dof_set;
DoFTools::extract_locally_active_dofs (dofh, dof_set);
TrilinosWrappers::MPI::Vector x_rel;
x_rel.reinit(relevant_set, MPI_COMM_WORLD);
- x_rel = 0;
x_rel.compress();
-
+
ConstraintMatrix cm(relevant_set);
DoFTools::make_hanging_node_constraints (dofh, cm);
cm.close ();
-
+
cm.distribute(x);
x_rel = x;
//x.print(std::cout);
-
+
// x_rel.print(std::cout);
TrilinosWrappers::Vector x_dub;
x_dub.reinit(dof_set.size());
-
+
x_dub = x_rel;
{
std::stringstream out;
out << "**** proc " << myid << std::endl;
x_dub.print (out);
-
+
if (myid==0)
deallog << out.str() << std::endl;
else
MPI_Send((void*)out.str().c_str(),out.str().size()+1, MPI_CHAR, 0, 1, MPI_COMM_WORLD);
}
-
+
if (myid == 0)
- {
+ {
for (unsigned int i = 1;i < numproc;++i)
{
MPI_Status status;
// $Id$
// Version: $Name$
//
-// Copyright (C) 2009, 2010, 2011 by the deal.II authors
+// Copyright (C) 2009, 2010, 2011, 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
trans.interpolate(x);
x_rel.reinit(relevant_set, MPI_COMM_WORLD);
- x_rel = 0;
x_rel.compress();
ConstraintMatrix cm(relevant_set);