]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merge remote-tracking branch 'mirror/trunk' into branch_cmake
authormaier <maier@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 Aug 2012 09:43:04 +0000 (09:43 +0000)
committermaier <maier@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 9 Aug 2012 09:43:04 +0000 (09:43 +0000)
Conflicts:
deal.II/source/Makefile

git-svn-id: https://svn.dealii.org/branches/branch_cmake@25798 0785d39b-7218-0410-832d-ea1e28bc413d

32 files changed:
deal.II/common/Make.global_options.in
deal.II/doc/news/changes.h
deal.II/doc/publications/index.html
deal.II/doc/readme.html
deal.II/examples/step-32/doc/intro.dox
deal.II/examples/step-42/doc/intro-step-42.tex [new file with mode: 0644]
deal.II/include/deal.II/integrators/laplace.h
deal.II/include/deal.II/lac/petsc_matrix_base.h
deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h
deal.II/include/deal.II/lac/petsc_sparse_matrix.h
deal.II/include/deal.II/lac/trilinos_vector.h
deal.II/include/deal.II/lac/trilinos_vector_base.h
deal.II/include/deal.II/numerics/point_value_history.h
deal.II/source/lac/petsc_matrix_base.cc
deal.II/source/lac/petsc_parallel_sparse_matrix.cc
deal.II/source/lac/petsc_sparse_matrix.cc
deal.II/source/numerics/dof_output_operator.cc [moved from deal.II/source/numerics/operator.cc with 96% similarity]
deal.II/source/numerics/dof_output_operator.inst.in [moved from deal.II/source/numerics/operator.inst.in with 100% similarity]
deal.II/source/numerics/point_value_history.cc
tests/deal.II/point_value_history_01.cc
tests/deal.II/point_value_history_01/cmp/generic
tests/deal.II/point_value_history_02.cc [new file with mode: 0644]
tests/deal.II/point_value_history_02/cmp/generic [new file with mode: 0644]
tests/deal.II/point_value_history_03.cc [new file with mode: 0644]
tests/deal.II/point_value_history_03/cmp/generic [new file with mode: 0644]
tests/integrators/laplacian_01.cc [new file with mode: 0644]
tests/integrators/laplacian_01/cmp/generic [new file with mode: 0644]
tests/mpi/constraint_matrix_trilinos_bug.cc
tests/mpi/p4est_2d_constraintmatrix_03.cc
tests/mpi/p4est_2d_constraintmatrix_04.cc
tests/mpi/p4est_3d_constraintmatrix_02.cc
tests/mpi/p4est_3d_constraintmatrix_03.cc

index 3a0db5f6b1af5b9440d011d42d9f4e083f6f9d40..80928d3ba8ae29ed3ca861f314682b6181fa561b 100644 (file)
@@ -157,7 +157,10 @@ endif
 # 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)
@@ -167,33 +170,27 @@ ifeq ($(USE_CONTRIB_PETSC),yes)
     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
 
index 0cbea7afb43fa220d42f441faaadea0a4ad3a549..9fbc5f42bc3bee726cffe05f486af5832f368f08 100644 (file)
@@ -24,6 +24,18 @@ inconvenience this causes.
 </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
@@ -50,6 +62,18 @@ used to store boundary indicators internally.
 
 
 <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
index 0cf3ba50812807244ea4684242d1c04b88bd955a..e10389e24a692a9bb4242fa8a905f7a814a83a22 100644 (file)
           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&#969 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&#969; 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>
index 2524fbc5428d46f2d6eccbe6af155fd6c8c83557..10c07e859aa56e8a721484fb9f7e6fd77ede76e1 100644 (file)
@@ -55,9 +55,9 @@
     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
index f26e7ef41ae958cb52c726823eeb7b15d8502b3d..143e8f6bfbe03ebd92ab11d271e294273d078ea3 100644 (file)
@@ -655,6 +655,10 @@ processor potentially has enough memory to hold the entries of the combined
 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
diff --git a/deal.II/examples/step-42/doc/intro-step-42.tex b/deal.II/examples/step-42/doc/intro-step-42.tex
new file mode 100644 (file)
index 0000000..a7e1fce
--- /dev/null
@@ -0,0 +1,182 @@
+\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
index bebeb553412420bb3d1a9d07137f4dd3a01c0ba0..802943abfa23bdbf33a177d5557799c9dd30a371 100644 (file)
@@ -69,6 +69,71 @@ namespace LocalIntegrators
             }
         }
     }
+
+/**
+ * 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[
@@ -110,9 +175,8 @@ namespace LocalIntegrators
     }
 
 /**
- * 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]
@@ -131,40 +195,38 @@ namespace LocalIntegrators
       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]
@@ -174,36 +236,42 @@ namespace LocalIntegrators
  * <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);
+                }
           }
       }
 
@@ -279,6 +347,136 @@ namespace LocalIntegrators
            }
        }
     }
+
+/**
+ * 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
@@ -291,6 +489,7 @@ namespace LocalIntegrators
  * <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
  */
index e35ecd55bfb26137d2b23477c3deb00d77b0b912..727eee8567714bf134be7d103c10e6c6bd9fa7ed 100644 (file)
@@ -952,68 +952,6 @@ namespace PETScWrappers
       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
index 6978daf8eda5b880820762bbd7fa12027366f5c8..dae469c5286125ee57e78772945bda72c537abbd 100644 (file)
@@ -18,6 +18,7 @@
 
 #  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
@@ -32,6 +33,8 @@ namespace PETScWrappers
   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
@@ -109,6 +112,7 @@ namespace PETScWrappers
     class SparseMatrix : public MatrixBase
     {
       public:
+
                                          /**
                                           * A structure that describes some of
                                           * the traits of this class in terms
@@ -381,6 +385,50 @@ namespace PETScWrappers
                         << "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:
 
                                          /**
index b3f21fe83c9821bd3e6854af8231d049b22922ee..aebee7087b0a7e6ed2a372036da01b4ecea4d3ba 100644 (file)
@@ -44,6 +44,7 @@ namespace PETScWrappers
   class SparseMatrix : public MatrixBase
   {
     public:
+
                                        /**
                                         * A structure that describes some of
                                         * the traits of this class in terms of
@@ -275,6 +276,49 @@ namespace PETScWrappers
                                         */
       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:
 
                                        /**
index ac6cb5582b847109e12c4ce9c66f26ca93f2ddff..f70303b23b2b9180ebd1e54ac802dc7a88815c21 100644 (file)
@@ -443,10 +443,13 @@ namespace TrilinosWrappers
                          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 &parallel_partitioning,
                      const MPI_Comm &communicator = MPI_COMM_WORLD,
@@ -658,19 +661,19 @@ namespace TrilinosWrappers
 
                                        /**
                                         * 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,
index c3d422bcfce8472e275fc29ee754018c413a49c4..549bc76afbd10f389f9995ea860aa6b81ae1925a 100644 (file)
@@ -1243,6 +1243,9 @@ namespace TrilinosWrappers
   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());
 
@@ -1260,6 +1263,10 @@ namespace TrilinosWrappers
   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()));
 
@@ -1273,6 +1280,10 @@ namespace TrilinosWrappers
   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()));
 
@@ -1616,6 +1627,9 @@ namespace TrilinosWrappers
   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();
@@ -1630,6 +1644,9 @@ namespace TrilinosWrappers
   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()));
 
@@ -1648,6 +1665,9 @@ namespace TrilinosWrappers
                    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(),
@@ -1668,6 +1688,9 @@ namespace TrilinosWrappers
   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()));
 
@@ -1686,6 +1709,9 @@ namespace TrilinosWrappers
                     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()));
 
@@ -1707,6 +1733,9 @@ namespace TrilinosWrappers
                     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(),
@@ -1733,6 +1762,9 @@ namespace TrilinosWrappers
                     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(),
@@ -1762,6 +1794,9 @@ namespace TrilinosWrappers
   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()));
 
@@ -1776,6 +1811,9 @@ namespace TrilinosWrappers
   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.
@@ -1804,6 +1842,9 @@ namespace TrilinosWrappers
                    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()));
 
index 09ed874c078854780ddb32ca5a395204d7e9b12e..b14c57586b0d8c29cee200028ee14e59b099a2f0 100644 (file)
@@ -143,14 +143,19 @@ namespace internal
  * 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
@@ -304,20 +309,21 @@ class PointValueHistory
                                      /**
                                       * 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
@@ -330,21 +336,21 @@ class PointValueHistory
     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);
 
 
@@ -353,7 +359,7 @@ class PointValueHistory
                                       * 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
@@ -372,7 +378,7 @@ class PointValueHistory
                                       * 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.
@@ -427,7 +433,7 @@ class PointValueHistory
                                       * 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
@@ -649,8 +655,8 @@ class PointValueHistory
                                       * 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
@@ -695,7 +701,7 @@ class PointValueHistory
     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.
@@ -737,12 +743,12 @@ class PointValueHistory
                                       */
     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
@@ -807,15 +813,15 @@ class PointValueHistory
     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 ();
 };
 
index d074a21ce8106e0cff3b96b0379b53a2bc749cd8..e0e3220c2f63bb435c1ff6eebc89cb58337e58dc 100644 (file)
@@ -549,28 +549,6 @@ namespace PETScWrappers
   }
 
 
-
-  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,
index ef7b6c3450497a683d405d9c0a05143fcc0772db..dddd2f579a7ac45e20200ebfd31905fbe70e892d 100644 (file)
@@ -613,8 +613,6 @@ namespace PETScWrappers
         }
     }
 
-
-
     // explicit instantiations
     //
     template
@@ -679,6 +677,25 @@ namespace PETScWrappers
                              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;
+    }
+    
   }
 }
 
index 6e15c4d5bb89663cb0e50edc897fecaa4a4d9e36..c69443649718ae4b7d16ee892eb1f31789513bf0 100644 (file)
@@ -25,6 +25,7 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace PETScWrappers
 {
+
   SparseMatrix::SparseMatrix ()
   {
     const int m=0, n=0, n_nonzero_per_row=0;
@@ -351,6 +352,23 @@ namespace PETScWrappers
   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;
+  }
 }
 
 
similarity index 96%
rename from deal.II/source/numerics/operator.cc
rename to deal.II/source/numerics/dof_output_operator.cc
index 213b2abb0bd4e4b367e6f8bb3b5d02f183d432e1..8fd7ff76ad750dfd85b07aa91eac83b37d75fea2 100644 (file)
@@ -31,7 +31,7 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace Algorithms
 {
-#include "operator.inst"
+#include "dof_output_operator.inst"
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 57fdb4fb287a2570ef931b6eb76867fb7a0b0201..a708e95cd4d694fa9e27d68753eb3a0b52c651e3 100644 (file)
@@ -24,6 +24,8 @@
 
 #include <deal.II/numerics/point_value_history.h>
 
+#include <algorithm>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -101,7 +103,7 @@ PointValueHistory<dim>::PointValueHistory (const PointValueHistory & point_value
   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;
 
@@ -133,7 +135,7 @@ PointValueHistory<dim>::operator= (const PointValueHistory & point_value_history
   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;
 
@@ -234,6 +236,14 @@ void PointValueHistory<dim>
 
                    // 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);
@@ -276,7 +286,7 @@ void PointValueHistory<dim>
                    // 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
@@ -298,8 +308,9 @@ void PointValueHistory<dim>
     {
                        // 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));
         }
@@ -372,6 +383,13 @@ void PointValueHistory<dim>
 
                    // 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);
@@ -414,8 +432,9 @@ void PointValueHistory<dim>
         {
                        // 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));
             }
@@ -430,7 +449,7 @@ void PointValueHistory<dim>
 
 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
@@ -440,16 +459,13 @@ void PointValueHistory<dim>
   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
@@ -461,7 +477,8 @@ void PointValueHistory<dim>
                    // 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;
@@ -473,13 +490,8 @@ template <int dim>
 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);
 }
 
 
@@ -491,9 +503,10 @@ void PointValueHistory<dim>
   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;
 }
@@ -562,27 +575,29 @@ void PointValueHistory<dim>
                            // 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++;
+        }
     }
   }
 }
@@ -731,18 +746,22 @@ void PointValueHistory<dim>
         {
             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
@@ -773,9 +792,9 @@ void PointValueHistory<dim>
   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
@@ -783,16 +802,12 @@ void PointValueHistory<dim>
                    // 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());
@@ -801,16 +816,17 @@ void PointValueHistory<dim>
                       // 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++;
+      }
     }
   }
 }
@@ -975,20 +991,23 @@ void PointValueHistory<dim>
 
           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 << "> ";
                     }
                 }
             }
@@ -1008,11 +1027,12 @@ void PointValueHistory<dim>
 
               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";
@@ -1183,25 +1203,29 @@ void PointValueHistory<dim>
 
   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
@@ -1285,17 +1309,17 @@ template <int dim>
 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;
 }
 
index d7c398ea818d5934be3326c50ec4ff49abb2f358..c4c241819c95a11524f16033f10edfe8881d9564 100644 (file)
@@ -83,8 +83,8 @@ void TestPointValueHistory<dim>::run()
     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);
 
@@ -187,9 +187,9 @@ void TestPointValueHistory<dim>::run()
 
         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;
@@ -219,39 +219,39 @@ void TestPointValueHistory<dim>::run()
 
     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;
+    }
 }
 
 
index 76e9c6c7e08627f6bd08935154afb23e170d755f..3698fc334d8841dc5d5ab8dd19ee142c35cb60f0 100644 (file)
@@ -3,6 +3,8 @@
 
 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)
@@ -12,7 +14,8 @@ Geometric Data
 
 
 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***
 
@@ -20,6 +23,8 @@ Solution: 3 : 1
 
 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)
@@ -59,8 +64,9 @@ Geometric Data
 
 
 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***
 
@@ -68,6 +74,8 @@ Solution: 18 : 10
 
 Closed: 1
 Cleared: 0
+Triangulation_changed: 0
+Have_dof_handler: 0
 Geometric Data
 No points stored currently
 
diff --git a/tests/deal.II/point_value_history_02.cc b/tests/deal.II/point_value_history_02.cc
new file mode 100644 (file)
index 0000000..272c171
--- /dev/null
@@ -0,0 +1,437 @@
+//----------------------------  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();
+}
diff --git a/tests/deal.II/point_value_history_02/cmp/generic b/tests/deal.II/point_value_history_02/cmp/generic
new file mode 100644 (file)
index 0000000..5555f1d
--- /dev/null
@@ -0,0 +1,271 @@
+
+***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::
diff --git a/tests/deal.II/point_value_history_03.cc b/tests/deal.II/point_value_history_03.cc
new file mode 100644 (file)
index 0000000..040dcda
--- /dev/null
@@ -0,0 +1,391 @@
+//----------------------------  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();
+}
diff --git a/tests/deal.II/point_value_history_03/cmp/generic b/tests/deal.II/point_value_history_03/cmp/generic
new file mode 100644 (file)
index 0000000..7ffc74d
--- /dev/null
@@ -0,0 +1,231 @@
+
+***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::
diff --git a/tests/integrators/laplacian_01.cc b/tests/integrators/laplacian_01.cc
new file mode 100644 (file)
index 0000000..10bbb9b
--- /dev/null
@@ -0,0 +1,247 @@
+//--------------------------------------------------------------------
+//    $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);
+  
+}
diff --git a/tests/integrators/laplacian_01/cmp/generic b/tests/integrators/laplacian_01/cmp/generic
new file mode 100644 (file)
index 0000000..059b0eb
--- /dev/null
@@ -0,0 +1,180 @@
+
+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
index 6a4761b0c9522e324b2a7ccffc10421024b3a06f..0d50c177f077080a9c28618951b112b1662efd1b 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -97,9 +97,9 @@ void test()
   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)
         {
@@ -121,7 +121,7 @@ void test()
   dofh.distribute_dofs (fe);
 
   IndexSet owned_set = dofh.locally_owned_dofs();
-  
+
   IndexSet dof_set;
   DoFTools::extract_locally_active_dofs (dofh, dof_set);
 
@@ -135,13 +135,12 @@ void test()
 
   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,
@@ -149,17 +148,17 @@ void test()
                                            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)
@@ -176,7 +175,7 @@ void test()
                                     //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;
 */
index 450fd20dd3daeea935bd9533c4cd338add38bf4d..d4cbf940b4931a262b818f298546f82dcfed58ac 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -147,7 +147,6 @@ void test()
 
   TrilinosWrappers::MPI::Vector x_rel;
   x_rel.reinit(relevant_set, MPI_COMM_WORLD);
-  x_rel = 0;
   x_rel.compress();
 
   ConstraintMatrix cm(relevant_set);
index c3c8ad97bf77e512eb573c5cdbe1f4f05025f33a..3b68e67243ca007e5b533c088fe51e680b5f365c 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -264,7 +264,6 @@ void test()
       trans.interpolate(x);
 
       x_rel.reinit(relevant_set, MPI_COMM_WORLD);
-      x_rel = 0;
       x_rel.compress();
 
       ConstraintMatrix cm(relevant_set);
index 9fc5a23aaaf1ef3ab7c6644eb89983d50a6d039b..058763664fee5f116297a9ed92d91d5b959d9525 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -61,14 +61,14 @@ void test()
 
          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);
 
@@ -82,38 +82,37 @@ void test()
 
   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;
index c2153bdf81460c1822111439cdae1988079a22e7..6c14055139aeb218fbe751f49c2b1c00ca2f629c 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -260,7 +260,6 @@ void test()
       trans.interpolate(x);
 
       x_rel.reinit(relevant_set, MPI_COMM_WORLD);
-      x_rel = 0;
       x_rel.compress();
 
       ConstraintMatrix cm(relevant_set);

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.