Remove the first_cell/is_first_cell/current_update_flags() mechanism from FiniteElement.
What this patch does is address the issue that
`FiniteElement::get_data()` (as was previously the case with Mappings
as well) sets a flag in the returned `InternalData` object that
indicates that we are visiting a cell for the first time. `FEValues`
later sets this flag back once we're done with the first cell.
The way this is used is that the various finite element
implementations query `current_update_flags` in
`FE::fill_fe_*_values()`, which equals `update_each|update_once` if we
are on the first cell, and `update_each` on later cells. That means
that on the first cell we update fields in the output object that we
need to initialize only once, e.g., the values of shape functions on
quadrature points, but we won't do this again on later cells.
I don't particularly like this approach. It is opaque and took even me
a long time to reconstruct when I saw it (see #1305 for a
discussion). My preference is to do things we only need to do once in
`FE::get_data()`, and do things we need to do every time in
`FE::fill_fe_*_values()`. This kind of change was already implemented
for mappings in the patches references from #1305, and the current
pull request goes into the same direction.
The current PR only does one half of the necessary changes in order to
keep the patch reasonably self contained and readable. In particular,
what it really does is simply *always* do everything, by storing
`update_once|update_each` in `FE::InternalDataBase::update_each`. In
other words, it is the equivalent to telling finite element
implementations that we are always on the first cell. This requires
more work than in the previous state, although copying shape values
every time is arguably of small expense compared to transforming
gradients.
I will of course fix this up in upcoming patches. It requires
significant shuffling around. I'll open an issue in a minute
discussing the details. I'd simply like to propose this here already
as an incremental step forward, without making patch review too
onerous. The patch is mostly mechanical once you understand what it
does.
Matthias Maier [Thu, 29 Oct 2015 03:47:29 +0000 (22:47 -0500)]
Bugfix: Provide missing operator+ and operator- variants
Unfortunately, it is not possible any more to exploit an implicit
conversion from SymmetricTensor to Tensor due to the heavily templated
operator variants for Tensor.
Thus, provide mixed operator types. (And clean up symmetric_tensor.h
later).
Break dependency of MappingQ from MappingQGeneric.
As discussed in #1732, MappingQ doesn't have the is-a property with
regard to MappingQGeneric. Nor does it have it with regard to its previous
base class, MappingQ1, and this was already fixed previously (#1429). Rather, what
it should be is a has-a relationship with regard to both the MappingQGeneric
(used on cells on the boundary) and MappingQ1 (used for interior cells).
This patch implements this. The part of the patch for MappingQ itself is
relatively straightforward: have it store pointers to both Q1 and Qp
mappings, and where it currently dispatches to the base class
MappingQGeneric, just dispatch to the Qp object instead. We can
then break the inheritance from MappingQ
to MappingQGeneric, at the cost of reimplementing a couple of
trivial functions of the Mapping interface that were previously
provided by the MappingQGeneric class.
The bigger problems appear in the MappingQEulerian and
MappingC1 classes which previously overloaded functions
declared in MappingQGeneric and that are now no longer
available to the MappingQGeneric object we now keep a
pointer to, rather than have as base class. This is
worked around by having both MappingQEulerian and
MappingC1 declare their own internal classes derived
from MappingQGeneric that provide the now missing
functions.
Alistair Bentley [Tue, 27 Oct 2015 15:14:07 +0000 (11:14 -0400)]
Implemented BDM_grad_grad in 3D
* Added second derivative terms to polynomials_bdm class in 3D
* Added tests and output to verify 2nd derivatives are returning
proper values
* I carefully checked the output file to confirm the correct output
is being produced
Changed size_type in PreconditionChebyshev. Added static_asserts that
compare matrix and vector types in precondtioners. Updated
linear_operator_08 test to account for inability to wrap
IterativeInverse in a linear_operator.
Rename the 'fast' parameter to 'omit_zeroing_entries'.
This is a more accurate representation of what the parameter achieves,
and should be easier to understand. The fact that this mode is also
'fast' is sort of a side effect of the semantics the parameter
really carries.
In the Table class, use omit_default_initialization.
Alistair Bentley [Tue, 20 Oct 2015 19:57:50 +0000 (15:57 -0400)]
Correction to BDM degree method
* The documentation says this function should return one less than
the highest polynomial degree. This is consistent with the
notation used in the literature on BDM polynomials and is what I
would expect the function to return.
* However, the function actually returns the degree of the highest
polynomial degree because it calls PolynomialSpace.degree() (note:
in this context, PolynomialSpace is the complete polynomial space
of degree k, next then we add elements from k+1 to complete the
BDM space) which specifically states it does not return the
Maxmial degree of the polynomial, but the number of polynomials
given to the constructor.
* If I understand correctly, since the number of polynomials given
to the constructor includes constants, this returns degree+1.
* To correct this, one needs to be subtracted from the degree of the
polynomial_space.
* A typo in the 3D bdm_grad functions has also been corrected
* Made a correction in the compute method to reflect the change in the
degree method.
* Corrected out in the test polynomials_bdm_01.output file
* Edited some names in the test bdm_1.cc
Timo Heister [Fri, 23 Oct 2015 18:42:07 +0000 (14:42 -0400)]
fix bdm tests
- turns out we need one more iteration in the CG solve
- fix test output (this never worked before?!)
- speed up by skipping highest degree in debug mode
- small difference in release/debug output (dropping small entries), so
adding alternative .output file
Ammended native deal.II preconditioners and solvers for use in LinearOperators.
Improved: Native deal.II preconditioners can now be directly wrapped as
LinearOperators. This allows them to be used in the construction of more complex
linear solvers, such as where the approximate inverse to a matrix is to be
defined as a preconditioner to said matrix.
New: Added a LinearOperator test for deal.II SparseMatrix solvers and
preconditioners.
Added a test to check that the various combinations of deal.II solvers and
preconditioners (for the SparseMatrix type) work with the linear_operator and
inverse_operator functions. By 'work' one means that they both compile and also
produce a sensible answer.
Note that the preconditioners are expected to act as per usual (i.e. in the
preconditioner slot in inverse_operator) as well as a standard matrix type in
linear_operator. This is necessary as they will be used in this context later
when a schur_matrix_operator and approximate_schur_matrix_operator is developed.
Fixed: Reinitialised (zeroed) range vectors in the construction of an
inverse_operator before the solve calls. This corrects a bug wherein the
starting point of the solver was garbage leading to potential divergence in the
linear solver.