From: maier Date: Thu, 9 Aug 2012 09:43:04 +0000 (+0000) Subject: Merge remote-tracking branch 'mirror/trunk' into branch_cmake X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c41f1e0ef0fa5c89b7598fa740a53d6c383222b2;p=dealii-svn.git Merge remote-tracking branch 'mirror/trunk' into branch_cmake Conflicts: deal.II/source/Makefile git-svn-id: https://svn.dealii.org/branches/branch_cmake@25798 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/common/Make.global_options.in b/deal.II/common/Make.global_options.in index 3a0db5f6b1..80928d3ba8 100644 --- a/deal.II/common/Make.global_options.in +++ b/deal.II/common/Make.global_options.in @@ -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 diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0cbea7afb4..9fbc5f42bc 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -24,6 +24,18 @@ inconvenience this causes.

    +
  1. 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. +
    +(Wolfgang Bangerth, 2012/08/06) +
  2. 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.
      +
    1. +New: officially added support for clang 3.1 or newer. +
      +(Timo Heister and Wolfgang Bangerth, 2012/08/07) + +
    2. +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. +
      +(Timo Heister, 2012/08/07) +
    3. Fixed: On some 64-bit systems, we build deal.II with the -m64 flag but forgot to build UMFPACK with this flag as well, leading to diff --git a/deal.II/doc/publications/index.html b/deal.II/doc/publications/index.html index 0cf3ba5081..e10389e24a 100644 --- a/deal.II/doc/publications/index.html +++ b/deal.II/doc/publications/index.html @@ -393,6 +393,32 @@ Preprint, University of Heidelberg, Germany, 2011.
    4. +
    5. + R. Hartmann +
      + Higher-order and adaptive discontinuous Galerkin + methods with shock-capturing applied to transonic turbulent + delta wing flow + +
      + 50th AIAA Aerospace Sciences Meeting, AIAA 2012-459, 2012. +
    6. + +
    7. + R. Hartmann +
      + Higher-order and adaptive Discontinuous Galerkin + methods applied to turbulent delta wing flow +
      + 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. +
    8. +
    9. T. Heister, G. Rapin
      Efficient augmented Lagrangian-type preconditioning @@ -542,6 +568,23 @@ pp. 3186-3197, 2012.
    10. +
    11. M. Wallraff, T. Leicht, M. Lange-Hegermann +
      + Numerical flux functions for Reynolds-averaged Navier-Stokes + and kω turbulence model computations with a line-preconditioned + p-multigrid discontinuous Galerkin solver +
      + Int. J. Numer. Meth. Fluids, published online, 2012. +
    12. + +
    13. M. Wallraff, T. Leicht, M. Lange-Hegermann +
      + Numerical flux functions for RANS-kω computations with a + line-implicit p-multigrid DG solver +
      + 50th AIAA Aerospace Sciences Meeting, AIAA 2012-458, 2012. +
    14. +
    15. J. Willems
      Robust Multilevel Methods for General Symmetric Positive Definite Operators @@ -732,30 +775,6 @@ SIAM J. of Numerical Analysis, vol. 49, pp. 917-944, 2011.
    16. -
    17. R. Hartmann -
      - Higher-order and adaptive Discontinuous Galerkin - methods applied to turbulent delta wing flow -
      - 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. -
    18. - -
    19. - R. Hartmann -
      - Higher-order and adaptive discontinuous Galerkin - methods with shock-capturing applied to transonic turbulent - delta wing flow - -
      - 50thAIAA Aerospace Sciences Meeting, 2012. Submitted. -
    20. -
    21. R. Hartmann, J. Held, T. Leicht
      @@ -766,18 +785,6 @@ J. Comput. Phys., vol. 230, pp. 4268-4284, 2011.
    22. -
    23. R. Hartmann, J. Held, T. Leicht, F. Prill -
      - Error estimation and adaptive mesh refinement for aerodynamic flows - -
      - 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. -
    24. -
    25. T. Heister
      A Massively Parallel Finite Element Framework with @@ -918,7 +925,7 @@ refinement for aerodynamic flows.
      - 49th49th AIAA Aerospace Sciences Meeting, 2011. + 49th AIAA Aerospace Sciences Meeting AIAA 2011-212, 2011.
    26. @@ -1384,6 +1391,21 @@ Aerospace Science and Technology, vol. 14, pp. 512-529, 2010.
    27. +
    28. R. Hartmann, J. Held, T. Leicht, F. Prill +
      + Error estimation and adaptive mesh refinement for aerodynamic flows + +
      + 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, +
      + vol. 113 of Notes on numerical fluid mechanics and + multidisciplinary design, pp. 339-353, Springer, 2010. +
    29. +
    30. R. Hartmann, P. Houston
      diff --git a/deal.II/doc/readme.html b/deal.II/doc/readme.html index 2524fbc542..10c07e859a 100644 --- a/deal.II/doc/readme.html +++ b/deal.II/doc/readme.html @@ -55,9 +55,9 @@ developed and tested on the following platforms:
      • 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
      • 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
      • Mac OS X 10.2 or later, with gcc 4.0 or later (Please see the deal.II Wiki for installation issues);
      • Sun Solaris 2.5, 2.6, 7, 8, and 9 on Sparc, with diff --git a/deal.II/examples/step-32/doc/intro.dox b/deal.II/examples/step-32/doc/intro.dox index f26e7ef41a..143e8f6bfb 100644 --- a/deal.II/examples/step-32/doc/intro.dox +++ b/deal.II/examples/step-32/doc/intro.dox @@ -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 index 0000000000..a7e1fce28d --- /dev/null +++ b/deal.II/examples/step-42/doc/intro-step-42.tex @@ -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 diff --git a/deal.II/include/deal.II/integrators/laplace.h b/deal.II/include/deal.II/integrators/laplace.h index bebeb55341..802943abfa 100644 --- a/deal.II/include/deal.II/integrators/laplace.h +++ b/deal.II/include/deal.II/integrators/laplace.h @@ -69,6 +69,71 @@ namespace LocalIntegrators } } } + +/** + * Laplacian residual operator in weak form + * + * \f[ + * -(\nabla u, \nabla v) + * \f] + */ + template + inline void + cell_residual ( + Vector& result, + const FEValuesBase& fe, + const std::vector >& 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 + inline void + cell_residual ( + Vector& result, + const FEValuesBase& fe, + const VectorSlice > > >& 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;kF 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 F - * the vector + * Weak boundary condition for the Laplace operator by Nitsche, scalar + * version, namely on the face F 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& result, const FEValuesBase& fe, - const VectorSlice > >& input, - const VectorSlice > > >& Dinput, - const VectorSlice > >& data, + const std::vector& input, + const std::vector >& Dinput, + const std::vector& 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& n = fe.normal_vector(k); for (unsigned i=0;iF the vector + * Weak boundary condition for the Laplace operator by Nitsche, vector + * valued version, namely on the face F + * 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 * Dinput, respectively. g is the inhomogeneous * boundary value in the argument data. $\gamma$ is the usual * penalty parameter. + * + * @ingroup Integrators + * @author Guido Kanschat + * @date 2008, 2009, 2010 */ template void nitsche_residual ( Vector& result, const FEValuesBase& fe, - const std::vector& input, - const std::vector >& Dinput, - const std::vector& data, + const VectorSlice > >& input, + const VectorSlice > > >& Dinput, + const VectorSlice > >& 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& n = fe.normal_vector(k); for (unsigned i=0;i + void + ip_residual( + Vector& result1, + Vector& result2, + const FEValuesBase& fe1, + const FEValuesBase& fe2, + const std::vector& input1, + const std::vector >& Dinput1, + const std::vector& input2, + const std::vector >& 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& n = fe1.normal_vector(k); + + for (unsigned i=0;i& 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 + void + ip_residual( + Vector& result1, + Vector& result2, + const FEValuesBase& fe1, + const FEValuesBase& fe2, + const VectorSlice > >& input1, + const VectorSlice > > >& Dinput1, + const VectorSlice > >& input2, + const VectorSlice > > >& 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& n = fe1.normal_vector(k); + + for (unsigned i=0;i& 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 * hi is the length of Zi * orthogonal to the current face. * + * @ingroup Integrators * @author Guido Kanschat * @date 2010 */ diff --git a/deal.II/include/deal.II/lac/petsc_matrix_base.h b/deal.II/include/deal.II/lac/petsc_matrix_base.h index e35ecd55bf..727eee8567 100644 --- a/deal.II/include/deal.II/lac/petsc_matrix_base.h +++ b/deal.II/include/deal.II/lac/petsc_matrix_base.h @@ -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 diff --git a/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h b/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h index 6978daf8ed..dae469c528 100644 --- a/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/petsc_parallel_sparse_matrix.h @@ -18,6 +18,7 @@ # include # include +# include # include 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: /** diff --git a/deal.II/include/deal.II/lac/petsc_sparse_matrix.h b/deal.II/include/deal.II/lac/petsc_sparse_matrix.h index b3f21fe83c..aebee7087b 100644 --- a/deal.II/include/deal.II/lac/petsc_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/petsc_sparse_matrix.h @@ -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: /** diff --git a/deal.II/include/deal.II/lac/trilinos_vector.h b/deal.II/include/deal.II/lac/trilinos_vector.h index ac6cb5582b..f70303b23b 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector.h +++ b/deal.II/include/deal.II/lac/trilinos_vector.h @@ -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 + * fast determines whether the + * vector should be filled with zero + * (false) or left untouched (true). */ void reinit (const IndexSet ¶llel_partitioning, const MPI_Comm &communicator = MPI_COMM_WORLD, @@ -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 - * fast 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 + * fast 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, diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h index c3d422bcfc..549bc76afb 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector_base.h +++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h @@ -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 &indices, const std::vector &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 &indices, const ::dealii::Vector &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())); diff --git a/deal.II/include/deal.II/numerics/point_value_history.h b/deal.II/include/deal.II/numerics/point_value_history.h index 09ed874c07..b14c57586b 100644 --- a/deal.II/include/deal.II/numerics/point_value_history.h +++ b/deal.II/include/deal.II/numerics/point_value_history.h @@ -143,14 +143,19 @@ namespace internal * course at the expense of getting only an approximate result.) *
    * - * When recording a new mnemonic name, the user must supply a std::vector - * 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@ 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 &components = std::vector ()); + const std::vector &component_mask = std::vector ()); /** * 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 &component_names); - - /** - * Provide optional names for the - * independent values. These names will - * be used instead of "Indep_...", if - * supplied. - */ + const std::vector &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 &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 > > data_store; + /** + * Saves a component mask + * for each mnemonic. + */ + std::map > component_mask; - /** - * Saves a vector listing components - * associated with a mnemonic. - */ - std::map > 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 (); }; diff --git a/deal.II/source/lac/petsc_matrix_base.cc b/deal.II/source/lac/petsc_matrix_base.cc index d074a21ce8..e0e3220c2f 100644 --- a/deal.II/source/lac/petsc_matrix_base.cc +++ b/deal.II/source/lac/petsc_matrix_base.cc @@ -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, diff --git a/deal.II/source/lac/petsc_parallel_sparse_matrix.cc b/deal.II/source/lac/petsc_parallel_sparse_matrix.cc index ef7b6c3450..dddd2f579a 100644 --- a/deal.II/source/lac/petsc_parallel_sparse_matrix.cc +++ b/deal.II/source/lac/petsc_parallel_sparse_matrix.cc @@ -613,8 +613,6 @@ namespace PETScWrappers } } - - // explicit instantiations // template @@ -679,6 +677,25 @@ namespace PETScWrappers const std::vector &, 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; + } + } } diff --git a/deal.II/source/lac/petsc_sparse_matrix.cc b/deal.II/source/lac/petsc_sparse_matrix.cc index 6e15c4d5bb..c694436497 100644 --- a/deal.II/source/lac/petsc_sparse_matrix.cc +++ b/deal.II/source/lac/petsc_sparse_matrix.cc @@ -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; + } } diff --git a/deal.II/source/numerics/operator.cc b/deal.II/source/numerics/dof_output_operator.cc similarity index 96% rename from deal.II/source/numerics/operator.cc rename to deal.II/source/numerics/dof_output_operator.cc index 213b2abb0b..8fd7ff76ad 100644 --- a/deal.II/source/numerics/operator.cc +++ b/deal.II/source/numerics/dof_output_operator.cc @@ -31,7 +31,7 @@ DEAL_II_NAMESPACE_OPEN namespace Algorithms { -#include "operator.inst" +#include "dof_output_operator.inst" } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/numerics/operator.inst.in b/deal.II/source/numerics/dof_output_operator.inst.in similarity index 100% rename from deal.II/source/numerics/operator.inst.in rename to deal.II/source/numerics/dof_output_operator.inst.in diff --git a/deal.II/source/numerics/point_value_history.cc b/deal.II/source/numerics/point_value_history.cc index 57fdb4fb28..a708e95cd4 100644 --- a/deal.II/source/numerics/point_value_history.cc +++ b/deal.II/source/numerics/point_value_history.cc @@ -24,6 +24,8 @@ #include +#include + DEAL_II_NAMESPACE_OPEN @@ -101,7 +103,7 @@ PointValueHistory::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::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 // 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 // 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 { // 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 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 (0)); } @@ -372,6 +383,13 @@ void PointValueHistory // 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 { // 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 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 (0)); } @@ -430,7 +449,7 @@ void PointValueHistory template void PointValueHistory -::add_field_name (const std::string &vector_name, const std::vector &components) +::add_field_name (const std::string &vector_name, const std::vector &mask) { // can't be closed to add additional points // or vectors @@ -440,16 +459,13 @@ void PointValueHistory AssertThrow (!triangulation_changed, ExcDoFHandlerChanged ()); - // Make a field_components list, if not supplied - std::vector temp_components = components; - if (temp_components.size() == 0) + // Make a component_mask, if not supplied + std::vector 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 (dof_handler->get_fe().n_components(), true); } - field_components.insert (std::pair > (vector_name, temp_components)); + component_mask.insert (std::pair > (vector_name, temp_mask)); // insert an empty vector of strings // to ensure each field has an entry @@ -461,7 +477,8 @@ void PointValueHistory // point_geometry_data.size() long std::pair > > 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 > vector_size (n_datastreams, std::vector (0)); pair_data.second = vector_size; @@ -473,13 +490,8 @@ template void PointValueHistory ::add_field_name(const std::string &vector_name, const unsigned int n_components) { - std::vector 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 temp_mask (n_components, true); + add_field_name (vector_name, temp_mask); } @@ -491,9 +503,10 @@ void PointValueHistory typename std::map >::iterator names = component_names_map.find(vector_name); Assert (names != component_names_map.end(), ExcMessage("vector_name not in class")); - typename std::map >::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 >::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 // to check vector_name is in the map. typename std::map > >::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 >::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 >::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 >::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 { typename std::map > >::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 >::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 >::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 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 // to check vector_name is in the map. typename std::map > >::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 >::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 >::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 >::iterator point = point_geometry_data.begin (); Vector value (dof_handler->get_fe().n_components()); @@ -801,16 +816,17 @@ void PointValueHistory // Make a Vector 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 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 >::iterator mask = component_mask.find(data_store_begin->first); + unsigned int n_stored = std::count(mask->second.begin(), mask->second.end(), true); std::vector 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 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 >::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 std::map > >::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 >::iterator components = field_components.find(vector_name); - Assert (components != field_components.end(), ExcMessage("vector_name not in class")); + typename std::map >::iterator mask = component_mask.find(vector_name); + Assert (mask != component_mask.end(), ExcMessage("vector_name not in class")); typename std::map >::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 void PointValueHistory ::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; } diff --git a/tests/deal.II/point_value_history_01.cc b/tests/deal.II/point_value_history_01.cc index d7c398ea81..c4c241819c 100644 --- a/tests/deal.II/point_value_history_01.cc +++ b/tests/deal.II/point_value_history_01.cc @@ -83,8 +83,8 @@ void TestPointValueHistory::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 finite_element(FE_Q (1 + 1), dim, FE_Q (1), 1); - DoFHandler dof_handler(triangulation); + // FESystem finite_element(FE_Q (1 + 1), dim, FE_Q (1), 1); // Already done! + // DoFHandler dof_handler(triangulation); dof_handler.distribute_dofs(finite_element); DoFRenumbering::Cuthill_McKee(dof_handler); @@ -187,9 +187,9 @@ void TestPointValueHistory::run() std::vector < std::vector > > selected_locations; node_monitor.get_points(selected_locations); - // write output to a file + // write output to a file Vector 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::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 + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include + +using namespace dealii; + + +template +class Postprocess : public DataPostprocessor +{ + 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 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 + // get_data_component_interpretation () const; +}; + +template +std::vector +Postprocess::get_names() const +{ + std::vector 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 +UpdateFlags +Postprocess::get_needed_update_flags () const +{ + return update_values | update_gradients | update_hessians; +} + +template +unsigned int +Postprocess::n_output_variables () const +{ + return 4; +} + + +template +void +Postprocess::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 +class TestPointValueHistory +{ +public: + TestPointValueHistory(); + void run(); + +private: + void output_results (unsigned int step, Vector solution) const; + + Triangulation triangulation; + FESystem finite_element; + DoFHandler dof_handler; + PointValueHistory test_copy; + std::vector > postprocessor_locations; +}; + + + + +template +TestPointValueHistory::TestPointValueHistory() : +finite_element(FE_Q (1 + 1), dim, FE_Q (1), 1), + +dof_handler(triangulation) +{ } + + + + +template +void TestPointValueHistory::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 finite_element(FE_Q (1 + 1), dim, FE_Q (1), 1); +// DoFHandler 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 block_component(dim + 1, 0); + block_component[dim] = 1; // component dim = pressure component! + DoFRenumbering::component_wise(dof_handler, block_component); + + // Vector + Vector 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 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 solution; +// solution.reinit(2); +// solution.block(0).reinit(n_u); +// solution.block(1).reinit(n_p); +// solution.collect_sizes(); +// +// BlockVector 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 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 quadrature_formula(finite_element.get_unit_support_points ()); + FEValues fe_values(finite_element, quadrature_formula, update_values | update_quadrature_points); // just need local_dof_indices and quadrature_points + + std::vector local_dof_indices(finite_element.dofs_per_cell); + std::vector > dof_locations(finite_element.dofs_per_cell); + + typename DoFHandler::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 node_monitor (dof_handler, n_inputs); + PointValueHistory 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 > > 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 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 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 component_mask (3, false); + component_mask[2] = true; + node_monitor.add_field_name("Pressure", component_mask); + component_mask = std::vector (3, false); + component_mask[1] = true; + node_monitor.add_field_name("Req_sol", component_mask); + component_mask = std::vector (4, true); + component_mask[3] = false; + node_monitor.add_field_name("Vector_out", component_mask); + component_mask = std::vector (4, false); + component_mask[3] = true; + node_monitor.add_field_name("Scalar_out", component_mask); + + std::vector 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_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 > > selected_locations; + node_monitor.get_support_locations(selected_locations); + Vector node_locations = node_monitor.mark_support_locations(); + QGauss 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 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 postprocessor; + QGauss postprocess_quadrature (2); + std::vector 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 +void TestPointValueHistory::output_results (unsigned int step, Vector solution) const +{ + std::vector solution_names (dim, "velocity"); + solution_names.push_back ("pressure"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, solution_names, + DataOut::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 index 0000000000..5555f1d1b6 --- /dev/null +++ b/tests/deal.II/point_value_history_02/cmp/generic @@ -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: +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 + +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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 index 0000000000..040dcda29f --- /dev/null +++ b/tests/deal.II/point_value_history_03.cc @@ -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 + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include + +using namespace dealii; + + +template +class Postprocess : public DataPostprocessor +{ + 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 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 + // get_data_component_interpretation () const; +}; + +template +std::vector +Postprocess::get_names() const +{ + std::vector names; + names.push_back ("X_gradient"); + names.push_back ("X_hessian"); + return names; +} + +template +UpdateFlags +Postprocess::get_needed_update_flags () const +{ + return update_values | update_gradients | update_hessians; +} + +template +unsigned int +Postprocess::n_output_variables () const +{ + return 2; +} + + + template + void + Postprocess::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 +class TestPointValueHistory +{ +public: + TestPointValueHistory(); + void run(); + +private: + void output_results (unsigned int step, Vector solution) const; + + Triangulation triangulation; + FE_Q finite_element; + DoFHandler dof_handler; + PointValueHistory test_copy; + std::vector > postprocessor_locations; +}; + + + + +template +TestPointValueHistory::TestPointValueHistory() : + finite_element(2), + dof_handler(triangulation) +{ } + + + + +template +void TestPointValueHistory::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 dof_handler(triangulation); + dof_handler.distribute_dofs(finite_element); + DoFRenumbering::Cuthill_McKee(dof_handler); + + // Vector + Vector 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 quadrature_formula(finite_element.get_unit_support_points ()); + FEValues fe_values(finite_element, quadrature_formula, update_values | update_quadrature_points); // just need local_dof_indices and quadrature_points + + std::vector local_dof_indices(finite_element.dofs_per_cell); + std::vector > dof_locations(finite_element.dofs_per_cell); + + typename DoFHandler::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 node_monitor (dof_handler, n_inputs); + PointValueHistory 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 > > 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 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 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 component_mask (1, true); + node_monitor.add_field_name("Req_sol", component_mask); + component_mask = std::vector (2, false); + component_mask[0] = true; + node_monitor.add_field_name("X_gradient", component_mask); + component_mask = std::vector (2, false); + component_mask[1] = true; + node_monitor.add_field_name("X_hessian", component_mask); + std::vector 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_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 > > selected_locations; + node_monitor.get_support_locations(selected_locations); + Vector node_locations = node_monitor.mark_support_locations(); + QGauss 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 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 postprocessor; + QGauss postprocess_quadrature (2); + std::vector 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 +void TestPointValueHistory::output_results (unsigned int step, Vector solution) const +{ + DataOut 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 index 0000000000..7ffc74d0f8 --- /dev/null +++ b/tests/deal.II/point_value_history_03/cmp/generic @@ -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: +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 + +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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::# +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 index 0000000000..10bbb9bb00 --- /dev/null +++ b/tests/integrators/laplacian_01.cc @@ -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 +#include + +#include +#include + +using namespace LocalIntegrators::Laplace; + +template +void test_cell(const FEValuesBase& fev) +{ + const unsigned int n = fev.dofs_per_cell; + unsigned int d=fev.get_fe().n_components(); + FullMatrix M(n,n); + cell_matrix(M,fev); + M.print(deallog,8); + + Vector u(n), v(n), w(n); + std::vector > > + ugrad(d,std::vector >(fev.n_quadrature_points)); + + std::vector indices(n); + for (unsigned int i=0;i +void test_boundary(const FEValuesBase& fev) +{ + const unsigned int n = fev.dofs_per_cell; + unsigned int d=fev.get_fe().n_components(); + FullMatrix M(n,n); + nitsche_matrix(M, fev, 17); + M.print(deallog,8); + + Vector u(n), v(n), w(n); + std::vector > + uval (d,std::vector(fev.n_quadrature_points)), + null_val(d,std::vector(fev.n_quadrature_points, 0.)); + std::vector > > + ugrad (d,std::vector >(fev.n_quadrature_points)); + + std::vector indices(n); + for (unsigned int i=0;i +void test_face(const FEValuesBase& fev1, + const FEValuesBase& 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 M11(n1,n1); + FullMatrix M12(n1,n2); + FullMatrix M21(n2,n1); + FullMatrix 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 u1(n1), v1(n1), w1(n1); + Vector u2(n2), v2(n2), w2(n2); + std::vector > + u1val (d,std::vector(fev1.n_quadrature_points)), + u2val (d,std::vector(fev2.n_quadrature_points)); + std::vector > > + u1grad(d,std::vector >(fev1.n_quadrature_points)), + u2grad(d,std::vector >(fev2.n_quadrature_points)); + + std::vector indices1(n1), indices2(n2); + for (unsigned int i=0;i +void +test_fe(Triangulation& tr, FiniteElement& fe) +{ + deallog << fe.get_name() << std::endl << "cell matrix" << std::endl; + QGauss quadrature(fe.tensor_degree()+1); + FEValues fev(fe, quadrature, update_gradients); + + typename Triangulation::cell_iterator cell1 = tr.begin(1); + fev.reinit(cell1); + test_cell(fev); + + QGauss face_quadrature(fe.tensor_degree()+1); + FEFaceValues fef1(fe, face_quadrature, update_values | update_gradients | update_normal_vectors); + for (unsigned int i=0;i::faces_per_cell;++i) + { + deallog << "boundary_matrix " << i << std::endl; + fef1.reinit(cell1, i); + test_boundary(fef1); + } + + FEFaceValues fef2(fe, face_quadrature, update_values | update_gradients); + typename Triangulation::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 +void +test(Triangulation& tr) +{ + FE_DGQ q1(1); + test_fe(tr, q1); + +// FE_DGQ q2(2); +// test_fe(tr, q2); + + FE_Nedelec 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 index 0000000000..059b0eb97e --- /dev/null +++ b/tests/integrators/laplacian_01/cmp/generic @@ -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 diff --git a/tests/mpi/constraint_matrix_trilinos_bug.cc b/tests/mpi/constraint_matrix_trilinos_bug.cc index 6a4761b0c9..0d50c177f0 100644 --- a/tests/mpi/constraint_matrix_trilinos_bug.cc +++ b/tests/mpi/constraint_matrix_trilinos_bug.cc @@ -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 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 velocity_mask (dim+1, true); - + velocity_mask[dim] = false; VectorTools::interpolate_boundary_values (dofh, @@ -149,17 +148,17 @@ void test() ZeroFunction(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; */ diff --git a/tests/mpi/p4est_2d_constraintmatrix_03.cc b/tests/mpi/p4est_2d_constraintmatrix_03.cc index 450fd20dd3..d4cbf940b4 100644 --- a/tests/mpi/p4est_2d_constraintmatrix_03.cc +++ b/tests/mpi/p4est_2d_constraintmatrix_03.cc @@ -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); diff --git a/tests/mpi/p4est_2d_constraintmatrix_04.cc b/tests/mpi/p4est_2d_constraintmatrix_04.cc index c3c8ad97bf..3b68e67243 100644 --- a/tests/mpi/p4est_2d_constraintmatrix_04.cc +++ b/tests/mpi/p4est_2d_constraintmatrix_04.cc @@ -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); diff --git a/tests/mpi/p4est_3d_constraintmatrix_02.cc b/tests/mpi/p4est_3d_constraintmatrix_02.cc index 9fc5a23aaa..058763664f 100644 --- a/tests/mpi/p4est_3d_constraintmatrix_02.cc +++ b/tests/mpi/p4est_3d_constraintmatrix_02.cc @@ -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 dofh(tr); static const FE_Q 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; diff --git a/tests/mpi/p4est_3d_constraintmatrix_03.cc b/tests/mpi/p4est_3d_constraintmatrix_03.cc index c2153bdf81..6c14055139 100644 --- a/tests/mpi/p4est_3d_constraintmatrix_03.cc +++ b/tests/mpi/p4est_3d_constraintmatrix_03.cc @@ -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);