From 40e0081292656589617aa47d2437f2b79b31da6c Mon Sep 17 00:00:00 2001 From: kronbichler Date: Wed, 10 Nov 2010 15:21:40 +0000 Subject: [PATCH] Change distributed and fe directories so that files are compiled only once and not for each dimension once. git-svn-id: https://svn.dealii.org/trunk@22676 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/fe/fe_abf.h | 8 - deal.II/include/deal.II/fe/fe_bdm.h | 5 +- deal.II/include/deal.II/fe/fe_dg_vector.h | 1 - deal.II/include/deal.II/fe/fe_nedelec.h | 6 +- .../include/deal.II/fe/fe_raviart_thomas.h | 2 - deal.II/include/deal.II/fe/fe_system.h | 88 +- deal.II/source/Makefile | 4 +- ...t.all_dimensions.cc => grid_refinement.cc} | 0 .../source/distributed/solution_transfer.cc | 32 +- .../distributed/solution_transfer.inst.in | 49 + deal.II/source/distributed/tria.cc | 196 ++-- deal.II/source/distributed/tria.inst.in | 64 ++ deal.II/source/fe/fe.cc | 207 ++--- deal.II/source/fe/fe.inst.in | 23 + deal.II/source/fe/fe_abf.cc | 142 +-- deal.II/source/fe/fe_abf.inst.in | 19 + deal.II/source/fe/fe_bdm.cc | 44 +- deal.II/source/fe/fe_bdm.inst.in | 19 + .../{fe_data.all_dimensions.cc => fe_data.cc} | 0 deal.II/source/fe/fe_dg_vector.cc | 5 +- deal.II/source/fe/fe_dg_vector.inst.in | 22 + deal.II/source/fe/fe_dgp.cc | 12 +- deal.II/source/fe/fe_dgp.inst.in | 23 + deal.II/source/fe/fe_dgp_monomial.cc | 49 +- deal.II/source/fe/fe_dgp_monomial.inst.in | 19 + deal.II/source/fe/fe_dgp_nonparametric.cc | 73 +- .../source/fe/fe_dgp_nonparametric.inst.in | 19 + deal.II/source/fe/fe_dgq.cc | 10 +- deal.II/source/fe/fe_dgq.inst.in | 25 + deal.II/source/fe/fe_face.cc | 18 +- deal.II/source/fe/fe_face.inst.in | 24 + deal.II/source/fe/fe_nedelec.cc | 438 +++++---- deal.II/source/fe/fe_nedelec.inst.in | 19 + deal.II/source/fe/fe_nothing.cc | 20 +- deal.II/source/fe/fe_nothing.inst.in | 19 + deal.II/source/fe/fe_poly.cc | 11 +- deal.II/source/fe/fe_poly.inst.in | 29 + deal.II/source/fe/fe_poly_tensor.cc | 213 +++-- deal.II/source/fe/fe_poly_tensor.inst.in | 22 + deal.II/source/fe/fe_q.cc | 34 +- deal.II/source/fe/fe_q.inst.in | 23 + deal.II/source/fe/fe_q_hierarchical.cc | 425 +++++---- deal.II/source/fe/fe_q_hierarchical.inst.in | 19 + deal.II/source/fe/fe_raviart_thomas.cc | 116 +-- deal.II/source/fe/fe_raviart_thomas.inst.in | 19 + deal.II/source/fe/fe_raviart_thomas_nodal.cc | 83 +- .../source/fe/fe_raviart_thomas_nodal.inst.in | 19 + deal.II/source/fe/fe_system.cc | 60 +- deal.II/source/fe/fe_system.inst.in | 22 + deal.II/source/fe/fe_tools.all_dimensions.cc | 302 ------ deal.II/source/fe/fe_tools.cc | 869 +++++++----------- deal.II/source/fe/fe_tools.inst.in | 522 +++++++++++ deal.II/source/fe/fe_values.cc | 104 +-- deal.II/source/fe/fe_values.inst.in | 53 +- deal.II/source/fe/mapping.cc | 10 +- deal.II/source/fe/mapping.inst.in | 23 + deal.II/source/fe/mapping_c1.cc | 23 +- deal.II/source/fe/mapping_c1.inst.in | 19 + deal.II/source/fe/mapping_cartesian.cc | 140 +-- deal.II/source/fe/mapping_cartesian.inst.in | 19 + deal.II/source/fe/mapping_q.cc | 44 +- deal.II/source/fe/mapping_q.inst.in | 23 + deal.II/source/fe/mapping_q1.cc | 121 ++- deal.II/source/fe/mapping_q1.inst.in | 25 + deal.II/source/fe/mapping_q1_eulerian.cc | 17 +- deal.II/source/fe/mapping_q1_eulerian.inst.in | 30 + deal.II/source/fe/mapping_q_eulerian.cc | 18 +- deal.II/source/fe/mapping_q_eulerian.inst.in | 33 + 68 files changed, 2815 insertions(+), 2379 deletions(-) rename deal.II/source/distributed/{grid_refinement.all_dimensions.cc => grid_refinement.cc} (100%) create mode 100644 deal.II/source/distributed/solution_transfer.inst.in create mode 100644 deal.II/source/distributed/tria.inst.in create mode 100644 deal.II/source/fe/fe.inst.in create mode 100644 deal.II/source/fe/fe_abf.inst.in create mode 100644 deal.II/source/fe/fe_bdm.inst.in rename deal.II/source/fe/{fe_data.all_dimensions.cc => fe_data.cc} (100%) create mode 100644 deal.II/source/fe/fe_dg_vector.inst.in create mode 100644 deal.II/source/fe/fe_dgp.inst.in create mode 100644 deal.II/source/fe/fe_dgp_monomial.inst.in create mode 100644 deal.II/source/fe/fe_dgp_nonparametric.inst.in create mode 100644 deal.II/source/fe/fe_dgq.inst.in create mode 100644 deal.II/source/fe/fe_face.inst.in create mode 100644 deal.II/source/fe/fe_nedelec.inst.in create mode 100644 deal.II/source/fe/fe_nothing.inst.in create mode 100644 deal.II/source/fe/fe_poly.inst.in create mode 100644 deal.II/source/fe/fe_poly_tensor.inst.in create mode 100644 deal.II/source/fe/fe_q.inst.in create mode 100644 deal.II/source/fe/fe_q_hierarchical.inst.in create mode 100644 deal.II/source/fe/fe_raviart_thomas.inst.in create mode 100644 deal.II/source/fe/fe_raviart_thomas_nodal.inst.in create mode 100644 deal.II/source/fe/fe_system.inst.in delete mode 100644 deal.II/source/fe/fe_tools.all_dimensions.cc create mode 100644 deal.II/source/fe/fe_tools.inst.in create mode 100644 deal.II/source/fe/mapping.inst.in create mode 100644 deal.II/source/fe/mapping_c1.inst.in create mode 100644 deal.II/source/fe/mapping_cartesian.inst.in create mode 100644 deal.II/source/fe/mapping_q.inst.in create mode 100644 deal.II/source/fe/mapping_q1.inst.in create mode 100644 deal.II/source/fe/mapping_q1_eulerian.inst.in create mode 100644 deal.II/source/fe/mapping_q_eulerian.inst.in diff --git a/deal.II/include/deal.II/fe/fe_abf.h b/deal.II/include/deal.II/fe/fe_abf.h index b2df64ce45..24c6fb817c 100644 --- a/deal.II/include/deal.II/fe/fe_abf.h +++ b/deal.II/include/deal.II/fe/fe_abf.h @@ -380,14 +380,6 @@ class FE_ABF : public FE_PolyTensor, dim> /*@}*/ -/* -------------- declaration of explicit specializations ------------- */ - -#ifndef DOXYGEN - -template <> -std::vector FE_ABF<1>::get_dpo_vector (const unsigned int); - -#endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/include/deal.II/fe/fe_bdm.h b/deal.II/include/deal.II/fe/fe_bdm.h index 8b91e325cb..8b6ebbc7ca 100644 --- a/deal.II/include/deal.II/fe/fe_bdm.h +++ b/deal.II/include/deal.II/fe/fe_bdm.h @@ -23,7 +23,8 @@ #include -using namespace dealii; + +DEAL_II_NAMESPACE_OPEN /** * The Brezzi-Douglas-Marini element. @@ -123,4 +124,6 @@ class FE_BDM std::vector > test_values; }; +DEAL_II_NAMESPACE_CLOSE + #endif diff --git a/deal.II/include/deal.II/fe/fe_dg_vector.h b/deal.II/include/deal.II/fe/fe_dg_vector.h index c351d7cff8..7ad6d58165 100644 --- a/deal.II/include/deal.II/fe/fe_dg_vector.h +++ b/deal.II/include/deal.II/fe/fe_dg_vector.h @@ -190,7 +190,6 @@ class FE_DGVector }; -/* -------------- declaration of explicit specializations ------------- */ DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/include/deal.II/fe/fe_nedelec.h b/deal.II/include/deal.II/fe/fe_nedelec.h index 79177cd035..065076d8ca 100644 --- a/deal.II/include/deal.II/fe/fe_nedelec.h +++ b/deal.II/include/deal.II/fe/fe_nedelec.h @@ -226,7 +226,7 @@ class FE_Nedelec : public FE_PolyTensor, dim> { const; virtual unsigned int memory_consumption () const; virtual FiniteElement * clone() const; - + private: /** * Only for internal use. Its @@ -333,7 +333,7 @@ class FE_Nedelec : public FE_PolyTensor, dim> { * integration. */ Table<2, double> boundary_weights; - + /** * Allow access from other * dimensions. @@ -345,8 +345,6 @@ class FE_Nedelec : public FE_PolyTensor, dim> { #ifndef DOXYGEN -template <> -std::vector FE_Nedelec<1>::get_dpo_vector (const unsigned int); template <> void FE_Nedelec<1>::initialize_restriction(); diff --git a/deal.II/include/deal.II/fe/fe_raviart_thomas.h b/deal.II/include/deal.II/fe/fe_raviart_thomas.h index 0b3c25321f..1e94f4949c 100644 --- a/deal.II/include/deal.II/fe/fe_raviart_thomas.h +++ b/deal.II/include/deal.II/fe/fe_raviart_thomas.h @@ -424,8 +424,6 @@ class FE_RaviartThomasNodal #ifndef DOXYGEN -template <> -std::vector FE_RaviartThomas<1>::get_dpo_vector (const unsigned int); template <> void FE_RaviartThomas<1>::initialize_restriction(); diff --git a/deal.II/include/deal.II/fe/fe_system.h b/deal.II/include/deal.II/fe/fe_system.h index 18c24781be..326e734277 100644 --- a/deal.II/include/deal.II/fe/fe_system.h +++ b/deal.II/include/deal.II/fe/fe_system.h @@ -143,7 +143,7 @@ DEAL_II_NAMESPACE_OPEN * see the documentation of Triangulation. * * @ingroup febase fe vector_valued - * + * * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2002, 2003, 2006, Ralf Hartmann 2001. */ template @@ -174,7 +174,7 @@ class FESystem : public FiniteElement FESystem (const FiniteElement &fe, const unsigned int n_elements); - /** + /** * Constructor for mixed * discretizations with two * base elements. @@ -184,7 +184,7 @@ class FESystem : public FiniteElement FESystem (const FiniteElement &fe1, const unsigned int n1, const FiniteElement &fe2, const unsigned int n2); - /** + /** * Constructor for mixed * discretizations with three * base elements. @@ -195,7 +195,7 @@ class FESystem : public FiniteElement const FiniteElement &fe2, const unsigned int n2, const FiniteElement &fe3, const unsigned int n3); - /** + /** * Constructor for mixed * discretizations with four * base elements. @@ -207,7 +207,7 @@ class FESystem : public FiniteElement const FiniteElement &fe3, const unsigned int n3, const FiniteElement &fe4, const unsigned int n4); - /** + /** * Constructor for mixed * discretizations with five * base elements. @@ -291,7 +291,7 @@ class FESystem : public FiniteElement virtual double shape_value_component (const unsigned int i, const Point &p, const unsigned int component) const; - + /** * Return the gradient of the * @p ith shape function at the @@ -398,7 +398,7 @@ class FESystem : public FiniteElement shape_grad_grad_component (const unsigned int i, const Point &p, const unsigned int component) const; - + /** * Return the matrix * interpolating from the given @@ -425,7 +425,7 @@ class FESystem : public FiniteElement get_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; - /** + /** * Number of different base * elements of this object. * @@ -458,7 +458,7 @@ class FESystem : public FiniteElement */ virtual const FiniteElement & base_element (const unsigned int index) const; - + /** * Check for non-zero values on a * face. @@ -484,7 +484,7 @@ class FESystem : public FiniteElement virtual Point unit_support_point (const unsigned int index) const; - + /** * Implementation of the * respective function in the @@ -495,10 +495,10 @@ class FESystem : public FiniteElement unit_face_support_point (const unsigned int index) const; /** - * @name Functions to support hp + * @name Functions to support hp * @{ */ - + /** * Return whether this element * implements its hanging node @@ -511,12 +511,12 @@ class FESystem : public FiniteElement * this function. */ virtual bool hp_constraints_are_implemented () const; - + /** * Return the matrix * interpolating from a face of * of one element to the face of - * the neighboring element. + * the neighboring element. * The size of the matrix is * then source.dofs_per_face times * this->dofs_per_face. @@ -537,13 +537,13 @@ class FESystem : public FiniteElement virtual void get_face_interpolation_matrix (const FiniteElement &source, FullMatrix &matrix) const; - + /** * Return the matrix * interpolating from a face of * of one element to the subface of - * the neighboring element. + * the neighboring element. * The size of the matrix is * then source.dofs_per_face times * this->dofs_per_face. @@ -625,7 +625,7 @@ class FESystem : public FiniteElement virtual std::vector > hp_quad_dof_identities (const FiniteElement &fe_other) const; - + /** * Return whether this element dominates * the one given as argument when they @@ -642,7 +642,7 @@ class FESystem : public FiniteElement FiniteElementDomination::Domination compare_for_face_domination (const FiniteElement &fe_other) const; //@} - + /** * Determine an estimate for the * memory consumption (in bytes) @@ -654,15 +654,15 @@ class FESystem : public FiniteElement * pointers to their base class, * rather than the class itself. */ - virtual unsigned int memory_consumption () const; - + virtual unsigned int memory_consumption () const; + protected: /** * Compute flags for initial * update only. */ virtual UpdateFlags update_once (const UpdateFlags flags) const; - + /** * Compute flags for update on * each cell. @@ -677,7 +677,7 @@ class FESystem : public FiniteElement * constructors of @p FESystem. */ virtual FiniteElement * clone() const; - + /** * Prepare internal data * structures and fill in values @@ -717,7 +717,7 @@ class FESystem : public FiniteElement * @p compute_fill that does the * work for all three * fill_fe*_values functions. - */ + */ virtual void fill_fe_face_values (const Mapping &mapping, const typename Triangulation::cell_iterator &cell, @@ -746,12 +746,12 @@ class FESystem : public FiniteElement typename Mapping::InternalDataBase &mapping_data, typename Mapping::InternalDataBase &fe_data, FEValuesData &data) const ; - + /** * Do the work for the three * fill_fe*_values functions. - * + * * Calls (among other things) * fill_fe_([sub]face)_values * of the base elements. Calls @@ -789,13 +789,13 @@ class FESystem : public FiniteElement * invalid. */ static const unsigned int invalid_face_number = numbers::invalid_unsigned_int; - + /** * Pairs of multiplicity and * element type. */ typedef std::pair *, unsigned int> ElementPair; - + /** * Pointer to underlying finite * element classes. @@ -836,7 +836,7 @@ class FESystem : public FiniteElement * constructor. */ void initialize_quad_dof_index_permutation (); - + /** * Helper function used in the constructor: * take a @p FiniteElementData object @@ -851,7 +851,7 @@ class FESystem : public FiniteElement static FiniteElementData multiply_dof_numbers (const FiniteElementData &fe_data, const unsigned int N); - + /** * Same as above for mixed elements * with two different sub-elements. @@ -918,7 +918,7 @@ class FESystem : public FiniteElement compute_restriction_is_additive_flags ( const FiniteElement &fe, const unsigned int N); - + /** * Same as above for mixed elements * with two different sub-elements. @@ -985,8 +985,8 @@ class FESystem : public FiniteElement compute_restriction_is_additive_flags ( const std::vector*> &fes, const std::vector &multiplicities); - - + + /** * Compute the non-zero vector * components of a composed @@ -1061,7 +1061,7 @@ class FESystem : public FiniteElement static std::vector > compute_nonzero_components (const std::vector*> &fes, const std::vector &multiplicities); - + /** * This function is simply * singled out of the @@ -1078,7 +1078,7 @@ class FESystem : public FiniteElement * Used by @p initialize. */ void build_cell_tables(); - + /** * Used by @p initialize. */ @@ -1099,7 +1099,7 @@ class FESystem : public FiniteElement template std::vector > hp_object_dof_identities (const FiniteElement &fe_other) const; - + /** * Usually: Fields of * cell-independent data. @@ -1122,7 +1122,7 @@ class FESystem : public FiniteElement * @p n_base_elements. */ InternalData (const unsigned int n_base_elements); - + /** * Destructor. Deletes all * @p InternalDatas whose @@ -1138,7 +1138,7 @@ class FESystem : public FiniteElement * be computed. */ bool compute_hessians; - + /** * Gives write-access to the * pointer to a @@ -1154,7 +1154,7 @@ class FESystem : public FiniteElement * pointer to a * @p InternalData of the * @p base_noth base element. - */ + */ typename FiniteElement::InternalDataBase & get_fe_data (const unsigned int base_no) const; @@ -1174,7 +1174,7 @@ class FESystem : public FiniteElement * pointer to a * @p FEValuesData for the * @p base_noth base element. - */ + */ FEValuesData & get_fe_values_data (const unsigned int base_no) const; /** @@ -1213,9 +1213,9 @@ class FESystem : public FiniteElement * of the sub-data objects. */ virtual void clear_first_cell (); - + private: - + /** * Pointers to * @p InternalData objects @@ -1263,12 +1263,6 @@ class FESystem : public FiniteElement }; -/* -------------- declaration of explicit specializations ------------- */ - - -template <> void FESystem<1>::initialize_unit_face_support_points (); -template <> void FESystem<1,2>::initialize_unit_face_support_points (); - DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/Makefile b/deal.II/source/Makefile index d07ffa33b5..7b93b2f61f 100644 --- a/deal.II/source/Makefile +++ b/deal.II/source/Makefile @@ -16,8 +16,8 @@ default all: debug optimized # list the dimension dependent and dimension independent directories -dim-dep-dirs := dofs fe hp multigrid numerics distributed -dim-indep-dirs := base lac grid +dim-dep-dirs := dofs hp multigrid numerics +dim-indep-dirs := base lac fe grid distributed all-dirs := $(dim-dep-dirs) $(dim-indep-dirs) # create lists of file names. diff --git a/deal.II/source/distributed/grid_refinement.all_dimensions.cc b/deal.II/source/distributed/grid_refinement.cc similarity index 100% rename from deal.II/source/distributed/grid_refinement.all_dimensions.cc rename to deal.II/source/distributed/grid_refinement.cc diff --git a/deal.II/source/distributed/solution_transfer.cc b/deal.II/source/distributed/solution_transfer.cc index 5ba92a7b85..4393bb6c5c 100644 --- a/deal.II/source/distributed/solution_transfer.cc +++ b/deal.II/source/distributed/solution_transfer.cc @@ -235,37 +235,7 @@ namespace parallel // explicit instantiations - -namespace parallel -{ - namespace distributed - { -#if deal_II_dimension > 1 - template class SolutionTransfer, DoFHandler >; - - -#ifdef DEAL_II_USE_PETSC - template class SolutionTransfer >; - - template class SolutionTransfer >; - - template class SolutionTransfer >; - - template class SolutionTransfer >; -#endif - -#ifdef DEAL_II_USE_TRILINOS - template class SolutionTransfer >; - - template class SolutionTransfer >; - - template class SolutionTransfer >; - - template class SolutionTransfer >; -#endif -#endif - } -} +#include "solution_transfer.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/distributed/solution_transfer.inst.in b/deal.II/source/distributed/solution_transfer.inst.in new file mode 100644 index 0000000000..f622eed773 --- /dev/null +++ b/deal.II/source/distributed/solution_transfer.inst.in @@ -0,0 +1,49 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + namespace parallel + \{ + namespace distributed + \{ +#if deal_II_dimension > 1 + template class SolutionTransfer, DoFHandler >; + + +#ifdef DEAL_II_USE_PETSC + template class SolutionTransfer >; + + template class SolutionTransfer >; + + template class SolutionTransfer >; + + template class SolutionTransfer >; +#endif + +#ifdef DEAL_II_USE_TRILINOS + template class SolutionTransfer >; + + template class SolutionTransfer >; + + template class SolutionTransfer >; + + template class SolutionTransfer >; +#endif +#endif + \} + \} + + } + diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index c90a8205f6..fb3b8c0628 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -193,7 +193,6 @@ namespace internal size_t (&connectivity_memory_used) (types<2>::connectivity * p4est); }; -#if deal_II_dimension == 2 int (&functions<2>::quadrant_compare) (const void *v1, const void *v2) = p4est_quadrant_compare; @@ -316,9 +315,7 @@ namespace internal size_t (&functions<2>::connectivity_memory_used) (types<2>::connectivity * p4est) = p4est_connectivity_memory_used; -#endif - template <> struct functions<3> { static @@ -444,8 +441,7 @@ namespace internal size_t (&connectivity_memory_used) (types<3>::connectivity * p4est); }; -#if deal_II_dimension == 3 - + int (&functions<3>::quadrant_compare) (const void *v1, const void *v2) = p8est_quadrant_compare; @@ -570,7 +566,6 @@ namespace internal size_t (&functions<3>::connectivity_memory_used) (types<3>::connectivity * p4est) = p8est_connectivity_memory_used; -#endif template @@ -1983,12 +1978,11 @@ namespace parallel -#if deal_II_dimension == 2 - - template + template <> void - Triangulation::copy_new_triangulation_to_p4est (internal::int2type<2>) + Triangulation<2,2>::copy_new_triangulation_to_p4est (internal::int2type<2>) { + const unsigned int dim = 2, spacedim = 2; Assert (this->n_cells(0) > 0, ExcInternalError()); Assert (this->n_levels() == 1, ExcInternalError()); @@ -2057,14 +2051,91 @@ namespace parallel /* user_pointer */ this); } + + + // TODO: This is a verbatim copy of the 2,2 + // case. However, we can't just specialize the + // dim template argument, but let spacedim open + template <> + void + Triangulation<2,3>::copy_new_triangulation_to_p4est (internal::int2type<2>) + { + const unsigned int dim = 2, spacedim = 3; + Assert (this->n_cells(0) > 0, ExcInternalError()); + Assert (this->n_levels() == 1, ExcInternalError()); + + // data structures that counts how many + // cells touch each vertex + // (vertex_touch_count), and which cells + // touch a given vertex (together with + // the local numbering of that vertex + // within the cells that touch it) + std::vector vertex_touch_count; + std::vector< + std::list< + std::pair::active_cell_iterator, + unsigned int> > > + vertex_to_cell; + get_vertex_to_cell_mappings (*this, + vertex_touch_count, + vertex_to_cell); + const internal::p4est::types<2>::locidx + num_vtt = std::accumulate (vertex_touch_count.begin(), + vertex_touch_count.end(), + 0); + + // now create a connectivity + // object with the right sizes + // for all arrays. set vertex + // information only in debug + // mode (saves a few bytes in + // optimized mode) + const bool set_vertex_info +#ifdef DEBUG + = true +#else + = false #endif + ; -#if deal_II_dimension == 3 + connectivity + = internal::p4est::functions<2>:: + connectivity_new ((set_vertex_info == true ? this->n_vertices() : 0), + this->n_cells(0), + this->n_vertices(), + num_vtt); - template + set_vertex_and_cell_info (*this, + vertex_touch_count, + vertex_to_cell, + coarse_cell_to_p4est_tree_permutation, + set_vertex_info, + connectivity); + + Assert (p4est_connectivity_is_valid (connectivity) == 1, + ExcInternalError()); + + // now create a forest out of the + // connectivity data structure + parallel_forest + = internal::p4est::functions<2>:: + new_forest (mpi_communicator, + connectivity, + /* minimum initial number of quadrants per tree */ 0, + /* minimum level of upfront refinement */ 0, + /* use uniform upfront refinement */ 1, + /* user_data_size = */ 0, + /* user_data_constructor = */ NULL, + /* user_pointer */ this); + } + + + + template <> void - Triangulation::copy_new_triangulation_to_p4est (internal::int2type<3>) + Triangulation<3,3>::copy_new_triangulation_to_p4est (internal::int2type<3>) { + const int dim = 3, spacedim = 3; Assert (this->n_cells(0) > 0, ExcInternalError()); Assert (this->n_levels() == 1, ExcInternalError()); @@ -2214,7 +2285,6 @@ namespace parallel /* user_pointer */ this); } -#endif template @@ -2686,6 +2756,7 @@ namespace parallel types::subdomain_id_t Triangulation::locally_owned_subdomain () const { + Assert (dim > 1, ExcNotImplemented()); return my_subdomain; } @@ -2802,7 +2873,7 @@ namespace parallel MPI_Comm Triangulation::get_communicator () const { - return mpi_communicator; + return mpi_communicator; } @@ -2902,48 +2973,75 @@ namespace parallel p4est_coarse_cell, attached_data_pack_callbacks); } + } + + // TODO: again problems with specialization in + // only one template argument + template <> + Triangulation<1,1>::Triangulation (MPI_Comm) + { + Assert (false, ExcNotImplemented()); + } + template <> + Triangulation<1,1>::~Triangulation () + { + Assert (false, ExcNotImplemented()); } + template <> + types::subdomain_id_t + Triangulation<1,1>::locally_owned_subdomain () const + { + Assert (false, ExcNotImplemented()); + return 0; + } -#if deal_II_dimension == 1 - template - Triangulation<1,spacedim>::Triangulation (MPI_Comm) + + template <> + MPI_Comm + Triangulation<1,1>::get_communicator () const + { + return MPI_COMM_WORLD; + } + + + + template <> + Triangulation<1,2>::Triangulation (MPI_Comm) { Assert (false, ExcNotImplemented()); } - template - Triangulation<1,spacedim>::~Triangulation () + template <> + Triangulation<1,2>::~Triangulation () { Assert (false, ExcNotImplemented()); } - template + template <> types::subdomain_id_t - Triangulation<1,spacedim>::locally_owned_subdomain () const + Triangulation<1,2>::locally_owned_subdomain () const { Assert (false, ExcNotImplemented()); return 0; } - template + template <> MPI_Comm - Triangulation<1,spacedim>::get_communicator () const + Triangulation<1,2>::get_communicator () const { return MPI_COMM_WORLD; } -#endif - } } @@ -2983,51 +3081,9 @@ namespace parallel -#ifdef DEAL_II_USE_P4EST - -// explicit instantiations -namespace internal -{ - namespace p4est - { -# if deal_II_dimension > 1 - template - void - init_quadrant_children - (const types::quadrant & p4est_cell, - types::quadrant (&p4est_children)[GeometryInfo::max_children_per_cell]); - - template - void - init_coarse_quadrant - (types::quadrant & quad); - - template - bool - quadrant_is_equal - (const types::quadrant & q1, - const types::quadrant & q2); - - template - bool - quadrant_is_ancestor - (const types::quadrant & q1, - const types::quadrant & q2); -#endif - } -} -namespace parallel -{ - namespace distributed - { - template class Triangulation; -# if deal_II_dimension < 3 - template class Triangulation; -# endif - } -} -#endif // DEAL_II_USE_P4EST +/*-------------- Explicit Instantiations -------------------------------*/ +#include "tria.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/distributed/tria.inst.in b/deal.II/source/distributed/tria.inst.in new file mode 100644 index 0000000000..f77f9f466b --- /dev/null +++ b/deal.II/source/distributed/tria.inst.in @@ -0,0 +1,64 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + +# ifdef DEAL_II_USE_P4EST + + namespace internal + \{ + namespace p4est + \{ +#if deal_II_dimension > 1 + template + void + init_quadrant_children + (const types::quadrant & p4est_cell, + types::quadrant (&p4est_children)[GeometryInfo::max_children_per_cell]); + + template + void + init_coarse_quadrant + (types::quadrant & quad); + + template + bool + quadrant_is_equal + (const types::quadrant & q1, + const types::quadrant & q2); + + template + bool + quadrant_is_ancestor + (const types::quadrant & q1, + const types::quadrant & q2); +# endif + \} + \} + namespace parallel + \{ + namespace distributed + \{ + template class Triangulation; +# if deal_II_dimension < 3 + template class Triangulation; +# endif + \} + \} + +# endif // DEAL_II_USE_P4EST + + } + diff --git a/deal.II/source/fe/fe.cc b/deal.II/source/fe/fe.cc index 32784952bc..b45135b670 100644 --- a/deal.II/source/fe/fe.cc +++ b/deal.II/source/fe/fe.cc @@ -33,11 +33,11 @@ DEAL_II_NAMESPACE_OPEN /*------------------------------- FiniteElement ----------------------*/ -template +template const double FiniteElement::fd_step_length = 1.0e-6; -template +template void FiniteElement:: InternalDataBase::initialize_2nd (const FiniteElement *element, @@ -59,7 +59,7 @@ InternalDataBase::initialize_2nd (const FiniteElement *element, // FEValues objects with slightly // shifted positions std::vector > diff_points (quadrature.size()); - + differences.resize(2*dim); for (unsigned int d=0; d *element, diff_points[i] = quadrature.point(i) - shift; const Quadrature minus_quad (diff_points, quadrature.get_weights()); differences[d+dim] = new FEValues (mapping, *element, - minus_quad, update_gradients); + minus_quad, update_gradients); } } -template +template FiniteElement::InternalDataBase::~InternalDataBase () { for (unsigned int i=0; i::InternalDataBase::~InternalDataBase () -template +template FiniteElement::FiniteElement ( const FiniteElementData &fe_data, const std::vector &r_i_a_f, @@ -119,7 +119,7 @@ FiniteElement::FiniteElement ( adjust_line_dof_index_for_line_orientation_table (dim == 3 ? this->dofs_per_line : 0), system_to_base_table(this->dofs_per_cell), - face_system_to_base_table(this->dofs_per_face), + face_system_to_base_table(this->dofs_per_face), component_to_base_table (this->components, std::make_pair(std::make_pair(0U, 0U), 0U)), restriction_is_additive_flags(r_i_a_f), @@ -141,7 +141,7 @@ FiniteElement::FiniteElement ( = const_cast&> (restriction_is_additive_flags); aux.resize(ndofs, restriction_is_additive_flags[0]); } - + if (nonzero_components.size() == 1 && ndofs > 1) { std::vector >& aux @@ -160,8 +160,8 @@ FiniteElement::FiniteElement ( std::bind2nd(std::not_equal_to(), 1U)) == n_nonzero_components_table.end()); - - + + Assert (restriction_is_additive_flags.size() == this->dofs_per_cell, ExcDimensionMismatch(restriction_is_additive_flags.size(), this->dofs_per_cell)); @@ -179,9 +179,9 @@ FiniteElement::FiniteElement ( Assert (n_nonzero_components_table[i] >= 1, ExcInternalError()); Assert (n_nonzero_components_table[i] <= this->n_components(), - ExcInternalError()); + ExcInternalError()); }; - + // initialize some tables in the // default way, i.e. if there is // only one (vector-)component; if @@ -194,12 +194,12 @@ FiniteElement::FiniteElement ( for (unsigned int j=0 ; jdofs_per_cell ; ++j) { system_to_component_table[j] = std::pair(0,j); - system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); + system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); } for (unsigned int j=0 ; jdofs_per_face ; ++j) { face_system_to_component_table[j] = std::pair(0,j); - face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); + face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j); } } // Fill with default value; may be @@ -229,14 +229,14 @@ FiniteElement::FiniteElement ( -template +template FiniteElement::~FiniteElement () {} -template +template double FiniteElement::shape_value (const unsigned int, const Point &) const @@ -247,7 +247,7 @@ FiniteElement::shape_value (const unsigned int, -template +template double FiniteElement::shape_value_component (const unsigned int, const Point &, @@ -259,7 +259,7 @@ FiniteElement::shape_value_component (const unsigned int, -template +template Tensor<1,dim> FiniteElement::shape_grad (const unsigned int, const Point &) const @@ -270,7 +270,7 @@ FiniteElement::shape_grad (const unsigned int, -template +template Tensor<1,dim> FiniteElement::shape_grad_component (const unsigned int, const Point &, @@ -282,7 +282,7 @@ FiniteElement::shape_grad_component (const unsigned int, -template +template Tensor<2,dim> FiniteElement::shape_grad_grad (const unsigned int, const Point &) const @@ -293,7 +293,7 @@ FiniteElement::shape_grad_grad (const unsigned int, -template +template Tensor<2,dim> FiniteElement::shape_grad_grad_component (const unsigned int, const Point &, @@ -304,7 +304,7 @@ FiniteElement::shape_grad_grad_component (const unsigned int, } -template +template void FiniteElement::reinit_restriction_and_prolongation_matrices ( const bool isotropic_restriction_only, @@ -314,7 +314,7 @@ FiniteElement::reinit_restriction_and_prolongation_matrices ( ref_case <= RefinementCase::isotropic_refinement; ++ref_case) { const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); - + for (unsigned int i=0; i::isotropic_refinement) @@ -350,7 +350,7 @@ FiniteElement::get_restriction_matrix (const unsigned int child, -template +template const FullMatrix & FiniteElement::get_prolongation_matrix (const unsigned int child, const RefinementCase &refinement_case) const @@ -372,7 +372,7 @@ FiniteElement::get_prolongation_matrix (const unsigned int child, //TODO:[GK] This is probably not the most efficient way of doing this. -template +template unsigned int FiniteElement::component_to_block_index (const unsigned int index) const { @@ -384,30 +384,21 @@ FiniteElement::component_to_block_index (const unsigned int index) } -template +template unsigned int FiniteElement::adjust_quad_dof_index_for_face_orientation (const unsigned int index, - const bool, - const bool, - const bool) const + const bool face_orientation, + const bool face_flip, + const bool face_rotation) const { // general template for 1D and 2D: not // implemented. in fact, the function // shouldn't even be called unless we are // in 3d, so throw an internal error - Assert (false, ExcInternalError()); - return index; -} - -#if deal_II_dimension == 3 + Assert (dim==3, ExcInternalError()); + if (dim < 3) + return index; -template <> -unsigned int -FiniteElement<3>::adjust_quad_dof_index_for_face_orientation (const unsigned int index, - const bool face_orientation, - const bool face_flip, - const bool face_rotation) const -{ // adjust dofs on 3d faces if the face is // flipped. note that we query a table that // derived elements need to have set up @@ -423,29 +414,20 @@ FiniteElement<3>::adjust_quad_dof_index_for_face_orientation (const unsigned int return index+adjust_quad_dof_index_for_face_orientation_table(index,4*face_orientation+2*face_flip+face_rotation); } -#endif - -template +template unsigned int FiniteElement::adjust_line_dof_index_for_line_orientation (const unsigned int index, - const bool) const + const bool line_orientation) const { // general template for 1D and 2D: do // nothing. Do not throw an Assertion, // however, in order to allow to call this // function in 2D as well - return index; -} - -#if deal_II_dimension == 3 + if (dim<3) + return index; -template <> -unsigned int -FiniteElement<3>::adjust_line_dof_index_for_line_orientation (const unsigned int index, - const bool line_orientation) const -{ Assert (indexdofs_per_line, ExcIndexRange(index,0,this->dofs_per_line)); Assert (adjust_line_dof_index_for_line_orientation_table.size()==this->dofs_per_line, ExcInternalError()); @@ -455,11 +437,9 @@ FiniteElement<3>::adjust_line_dof_index_for_line_orientation (const unsigned int return index+adjust_line_dof_index_for_line_orientation_table[index]; } -#endif - -template +template bool FiniteElement::prolongation_is_implemented () const { @@ -512,7 +492,7 @@ bool FiniteElement::isotropic_prolongation_is_implemented () const { const RefinementCase ref_case=RefinementCase::isotropic_refinement; - + for (unsigned int c=0; c::n_children(RefinementCase(ref_case)); ++c) { @@ -531,7 +511,7 @@ FiniteElement::isotropic_prolongation_is_implemented () const -template +template bool FiniteElement::isotropic_restriction_is_implemented () const { @@ -555,7 +535,7 @@ FiniteElement::isotropic_restriction_is_implemented () const -template +template bool FiniteElement::constraints_are_implemented (const internal::SubfaceCase &subface_case) const { @@ -567,7 +547,7 @@ FiniteElement::constraints_are_implemented (const internal::Subfac -template +template bool FiniteElement::hp_constraints_are_implemented () const { @@ -576,27 +556,27 @@ FiniteElement::hp_constraints_are_implemented () const -template +template const FullMatrix & FiniteElement::constraints (const internal::SubfaceCase &subface_case) const { Assert (subface_case==internal::SubfaceCase::case_isotropic, ExcConstraintsVoid()); Assert ((this->dofs_per_face == 0) || (interface_constraints.m() != 0), ExcConstraintsVoid()); - + if (dim==1) Assert ((interface_constraints.m()==0) && (interface_constraints.n()==0), ExcWrongInterfaceMatrixSize(interface_constraints.m(), interface_constraints.n())); - + return interface_constraints; } -template +template TableIndices<2> -FiniteElement::interface_constraints_size () const +FiniteElement::interface_constraints_size () const { switch (dim) { @@ -620,7 +600,7 @@ FiniteElement::interface_constraints_size () const -template +template void FiniteElement:: get_interpolation_matrix (const FiniteElement &, @@ -637,7 +617,7 @@ get_interpolation_matrix (const FiniteElement &, -template +template void FiniteElement:: get_face_interpolation_matrix (const FiniteElement &, @@ -651,10 +631,10 @@ get_face_interpolation_matrix (const FiniteElement &, typename FEE:: ExcInterpolationNotImplemented()); } - - -template + + +template void FiniteElement:: get_subface_interpolation_matrix (const FiniteElement &, @@ -668,10 +648,10 @@ get_subface_interpolation_matrix (const FiniteElement &, AssertThrow (false, typename FEE::ExcInterpolationNotImplemented()); } - -template + +template std::vector > FiniteElement:: hp_vertex_dof_identities (const FiniteElement &) const @@ -682,7 +662,7 @@ hp_vertex_dof_identities (const FiniteElement &) const -template +template std::vector > FiniteElement:: hp_line_dof_identities (const FiniteElement &) const @@ -693,7 +673,7 @@ hp_line_dof_identities (const FiniteElement &) const -template +template std::vector > FiniteElement:: hp_quad_dof_identities (const FiniteElement &) const @@ -704,7 +684,7 @@ hp_quad_dof_identities (const FiniteElement &) const -template +template FiniteElementDomination::Domination FiniteElement:: compare_for_face_domination (const FiniteElement &) const @@ -715,7 +695,7 @@ compare_for_face_domination (const FiniteElement &) const -template +template bool FiniteElement::operator == (const FiniteElement &f) const { @@ -726,7 +706,7 @@ FiniteElement::operator == (const FiniteElement &f) -template +template const std::vector > & FiniteElement::get_unit_support_points () const { @@ -742,7 +722,7 @@ FiniteElement::get_unit_support_points () const -template +template bool FiniteElement::has_support_points () const { @@ -751,7 +731,7 @@ FiniteElement::has_support_points () const -template +template const std::vector > & FiniteElement::get_generalized_support_points () const { @@ -766,7 +746,7 @@ FiniteElement::get_generalized_support_points () const -template +template bool FiniteElement::has_generalized_support_points () const { @@ -775,7 +755,7 @@ FiniteElement::has_generalized_support_points () const -template +template Point FiniteElement::unit_support_point (const unsigned index) const { @@ -788,7 +768,7 @@ FiniteElement::unit_support_point (const unsigned index) const -template +template const std::vector > & FiniteElement::get_unit_face_support_points () const { @@ -804,7 +784,7 @@ FiniteElement::get_unit_face_support_points () const -template +template bool FiniteElement::has_face_support_points () const { @@ -813,7 +793,7 @@ FiniteElement::has_face_support_points () const -template +template const std::vector > & FiniteElement::get_generalized_face_support_points () const { @@ -828,7 +808,7 @@ FiniteElement::get_generalized_face_support_points () const -template +template bool FiniteElement::has_generalized_face_support_points () const { @@ -837,7 +817,7 @@ FiniteElement::has_generalized_face_support_points () const -template +template Point FiniteElement::unit_face_support_point (const unsigned index) const { @@ -849,7 +829,7 @@ FiniteElement::unit_face_support_point (const unsigned index) cons } -template +template bool FiniteElement::has_support_on_face ( const unsigned int, @@ -859,7 +839,7 @@ FiniteElement::has_support_on_face ( } -template +template void FiniteElement::interpolate( std::vector& local_dofs, @@ -872,14 +852,14 @@ FiniteElement::interpolate( ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); Assert (this->n_components() == 1, ExcDimensionMismatch(this->n_components(), 1)); - + std::copy(values.begin(), values.end(), local_dofs.begin()); } -template +template void FiniteElement::interpolate( std::vector& local_dofs, @@ -893,7 +873,7 @@ FiniteElement::interpolate( ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); Assert (values[0].size() >= offset+this->n_components(), ExcDimensionMismatch(values[0].size(),offset+this->n_components())); - + for (unsigned int i=0;idofs_per_cell;++i) { const std::pair index @@ -905,7 +885,7 @@ FiniteElement::interpolate( -template +template void FiniteElement::interpolate( std::vector& local_dofs, @@ -918,7 +898,7 @@ FiniteElement::interpolate( ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); Assert (values.size() == this->n_components(), ExcDimensionMismatch(values.size(), this->n_components())); - + for (unsigned int i=0;idofs_per_cell;++i) { const std::pair index @@ -929,8 +909,7 @@ FiniteElement::interpolate( - -template +template unsigned int FiniteElement::memory_consumption () const { @@ -941,7 +920,7 @@ FiniteElement::memory_consumption () const MemoryConsumption::memory_consumption (system_to_component_table) + MemoryConsumption::memory_consumption (face_system_to_component_table) + MemoryConsumption::memory_consumption (system_to_base_table) + - MemoryConsumption::memory_consumption (face_system_to_base_table) + + MemoryConsumption::memory_consumption (face_system_to_base_table) + MemoryConsumption::memory_consumption (component_to_base_table) + MemoryConsumption::memory_consumption (restriction_is_additive_flags) + MemoryConsumption::memory_consumption (nonzero_components) + @@ -949,8 +928,6 @@ FiniteElement::memory_consumption () const } -#if deal_II_dimension == 1 - template<> void FiniteElement<1,2>::compute_2nd ( @@ -965,9 +942,7 @@ FiniteElement<1,2>::compute_2nd ( Assert(false, ExcNotImplemented()); } -#endif -#if deal_II_dimension == 2 template<> void @@ -983,9 +958,9 @@ FiniteElement<2,3>::compute_2nd ( Assert(false, ExcNotImplemented()); } -#endif -template + +template void FiniteElement::compute_2nd ( const Mapping &mapping, @@ -999,7 +974,7 @@ FiniteElement::compute_2nd ( & update_hessians, ExcInternalError()); -// make sure we have as many entries as there are nonzero components +// make sure we have as many entries as there are nonzero components // Assert (data.shape_hessians.size() == // std::accumulate (n_nonzero_components_table.begin(), // n_nonzero_components_table.end(), @@ -1098,7 +1073,7 @@ FiniteElement::compute_2nd ( left = fe_internal.differences[d1+dim] ->shape_grad_component(shape_index, q+offset, component); }; - + // compute the second // derivative from a // symmetric difference @@ -1106,7 +1081,7 @@ FiniteElement::compute_2nd ( for (unsigned int d=0; d::compute_2nd ( ExcInternalError()); mapping.transform (diff_quot[d], diff_quot2, mapping_internal, mapping_covariant); - + for (unsigned int q=0; q::compute_2nd ( -template +template std::vector FiniteElement::compute_n_nonzero_components ( const std::vector > &nonzero_components) @@ -1147,7 +1122,7 @@ FiniteElement::compute_n_nonzero_components ( /*------------------------------- FiniteElement ----------------------*/ -template +template typename Mapping::InternalDataBase * FiniteElement::get_face_data (const UpdateFlags flags, const Mapping &mapping, @@ -1159,7 +1134,7 @@ FiniteElement::get_face_data (const UpdateFlags flags, -template +template typename Mapping::InternalDataBase * FiniteElement::get_subface_data (const UpdateFlags flags, const Mapping &mapping, @@ -1171,7 +1146,7 @@ FiniteElement::get_subface_data (const UpdateFlags flags, -template +template const FiniteElement& FiniteElement::base_element(const unsigned index) const { @@ -1179,16 +1154,10 @@ FiniteElement::base_element(const unsigned index) const return *this; } -/*------------------------------- Explicit Instantiations -------------*/ - -template class FiniteElement; - - -#if deal_II_dimension != 3 -template class FiniteElement; -#endif +/*------------------------------- Explicit Instantiations -------------*/ +#include "fe.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe.inst.in b/deal.II/source/fe/fe.inst.in new file mode 100644 index 0000000000..656bd44e58 --- /dev/null +++ b/deal.II/source/fe/fe.inst.in @@ -0,0 +1,23 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FiniteElement; + +#if deal_II_dimension != 3 + template class FiniteElement; +#endif + } + diff --git a/deal.II/source/fe/fe_abf.cc b/deal.II/source/fe/fe_abf.cc index ea2c3a1ca9..6c14d14e8f 100644 --- a/deal.II/source/fe/fe_abf.cc +++ b/deal.II/source/fe/fe_abf.cc @@ -55,7 +55,7 @@ FE_ABF::FE_ABF (const unsigned int deg) { Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - + this->mapping_type = mapping_piola; // First, initialize the // generalized support points and @@ -74,7 +74,7 @@ FE_ABF::FE_ABF (const unsigned int deg) // From now on, the shape functions // will be the correct ones, not // the raw shape functions anymore. - + // Reinit the vectors of // restriction and prolongation // matrices to the right sizes. @@ -123,7 +123,7 @@ FE_ABF::get_name () const #else std::ostrstream namebuf; #endif - + namebuf << "FE_ABF<" << dim << ">(" << rt_order << ")"; #ifndef HAVE_STD_STRINGSTREAM @@ -147,38 +147,6 @@ FE_ABF::clone() const //--------------------------------------------------------------------------- -#if deal_II_dimension == 1 - -template -void -FE_ABF::initialize_support_points (const unsigned int deg) -{ - return; - - Assert (false, ExcNotImplemented()); - - QGauss cell_quadrature(deg+1); - const unsigned int n_interior_points - = (deg>0) ? cell_quadrature.size() : 0; - - this->generalized_support_points.resize (2 + n_interior_points); - - // Number of the point being entered - unsigned int current = 0; - - - if (deg==0) return; - - interior_weights.reinit(TableIndices<3>(2+n_interior_points, 0, dim)); - - for (unsigned int k=0;kgeneralized_support_points[current++] = cell_quadrature.point(k); - - Assert (current == this->generalized_support_points.size(), - ExcInternalError()); -} - -#else // Version for 2d and higher. See above for 1d version template @@ -219,14 +187,14 @@ FE_ABF::initialize_support_points (const unsigned int deg) if (dim>1) { QGauss face_points (deg+1); - TensorProductPolynomials legendre - = Polynomials::Legendre::generate_complete_basis(deg); + TensorProductPolynomials legendre = + Polynomials::Legendre::generate_complete_basis(deg); boundary_weights.reinit(n_face_points, legendre.n()); - + // Assert (face_points.size() == this->dofs_per_face, // ExcInternalError()); - + for (unsigned int k=0;kgeneralized_face_support_points[k] = face_points.point(k); @@ -266,13 +234,13 @@ FE_ABF::initialize_support_points (const unsigned int deg) } } } - + // Create Legendre basis for the // space D_xi Q_k if (deg>0) { std::vector* > polynomials(dim); - + for (unsigned int dd=0;dd > > poly(dim); @@ -282,9 +250,9 @@ FE_ABF::initialize_support_points (const unsigned int deg) polynomials[dd] = new AnisotropicPolynomials(poly); } - + interior_weights.reinit(TableIndices<3>(n_interior_points, polynomials[0]->n(), dim)); - + for (unsigned int k=0;kn();++i) @@ -292,13 +260,13 @@ FE_ABF::initialize_support_points (const unsigned int deg) interior_weights(k,i,d) = cell_quadrature.weight(k) * polynomials[d]->compute_value(i,cell_quadrature.point(k)); } - + for (unsigned int d=0;dgeneralized_support_points[current++] = cell_quadrature.point(k); @@ -307,7 +275,7 @@ FE_ABF::initialize_support_points (const unsigned int deg) // TODO: Here the canonical extension of the principle // behind the ABF elements is implemented. It is unclear, // if this really leads to the ABF spaces in 3D! - interior_weights_abf.reinit(TableIndices<3>(cell_quadrature.size(), + interior_weights_abf.reinit(TableIndices<3>(cell_quadrature.size(), polynomials_abf[0]->n() * dim, dim)); Tensor<1, dim> poly_grad; @@ -330,21 +298,7 @@ FE_ABF::initialize_support_points (const unsigned int deg) ExcInternalError()); } -#endif - -#if deal_II_dimension == 1 - -template -void -FE_ABF::initialize_restriction() -{ - unsigned int iso=RefinementCase::isotropic_refinement-1; - for (unsigned int i=0;i::max_children_per_cell;++i) - this->restriction[iso][i].reinit(0,0); -} - -#else // This function is the same Raviart-Thomas interpolation performed by // interpolate. Still, we cannot use interpolate, since it was written @@ -360,6 +314,13 @@ template void FE_ABF::initialize_restriction() { + if (dim==1) + { + unsigned int iso=RefinementCase::isotropic_refinement-1; + for (unsigned int i=0;i::max_children_per_cell;++i) + this->restriction[iso][i].reinit(0,0); + return; + } unsigned int iso=RefinementCase::isotropic_refinement-1; QGauss q_base (rt_order+1); const unsigned int n_face_points = q_base.size(); @@ -394,7 +355,7 @@ FE_ABF::initialize_restriction() const unsigned int child = GeometryInfo::child_cell_on_face( RefinementCase::isotropic_refinement, face, sub); - + // On a certain face, we must // compute the moments of ALL // fine level functions with @@ -426,9 +387,9 @@ FE_ABF::initialize_restriction() } } } - + if (rt_order==0) return; - + // Create Legendre basis for the // space D_xi Q_k. Here, we cannot // use the shape functions @@ -442,7 +403,7 @@ FE_ABF::initialize_restriction() polynomials[dd] = new AnisotropicPolynomials(poly); } - + QGauss q_cell(rt_order+1); const unsigned int start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; @@ -455,11 +416,11 @@ FE_ABF::initialize_restriction() for (unsigned int i = 0; i < this->dofs_per_cell; ++i) for (unsigned int d=0;dshape_value_component(i, q_cell.point(k), d); - + for (unsigned int child=0;child::max_children_per_cell;++child) { Quadrature q_sub = QProjector::project_to_child(q_cell, child); - + for (unsigned int k=0;kdofs_per_cell; ++i_child) for (unsigned int d=0;d::initialize_restriction() * polynomials[d]->compute_value(i_weight, q_sub.point(k)); } } - + for (unsigned int d=0;d -std::vector -FE_ABF<1>::get_dpo_vector (const unsigned int) -{ - Assert (false, ExcImpossibleInDim(1)); - return std::vector(); -} - -#endif template std::vector FE_ABF::get_dpo_vector (const unsigned int rt_order) { + if (dim == 1) + { + Assert (false, ExcImpossibleInDim(1)); + return std::vector(); + } + // the element is face-based (not // to be confused with George // W. Bush's Faith Based @@ -508,11 +462,11 @@ FE_ABF::get_dpo_vector (const unsigned int rt_order) // and then there are interior dofs const unsigned int interior_dofs = dim*(rt_order+1)*dofs_per_face; - + std::vector dpo(dim+1); dpo[dim-1] = dofs_per_face; dpo[dim] = interior_dofs; - + return dpo; } @@ -542,7 +496,7 @@ FE_ABF::update_each (const UpdateFlags flags) const | update_piola | update_JxW_values; if (flags & update_gradients) - out |= update_gradients | update_covariant_transformation + out |= update_gradients | update_covariant_transformation | update_piola | update_JxW_values; //TODO: Set update flags appropriately and figure out, how the second @@ -617,16 +571,16 @@ FE_ABF::has_support_on_face (const unsigned int shape_index, // these in a table return (face_index != GeometryInfo::opposite_face[shape_index]); } - + default: return true; }; }; - + default: // other rt_order return true; }; - + return true; } @@ -667,10 +621,10 @@ FE_ABF::interpolate( local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i) * values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]+offset); } - + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; - + for (unsigned int k=0;k::interpolate( // TODO: Check what the face_orientation, face_flip and face_rotation have to be in 3D unsigned int k = QProjector::DataSetDescriptor::face (face, false, false, false, n_face_points); for (unsigned int i=0; i::unit_normal_direction[face]+offset); } } @@ -727,7 +681,7 @@ FE_ABF::interpolate( local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i) * values[GeometryInfo::unit_normal_direction[face]][face*n_face_points+k]; } - + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; @@ -753,7 +707,7 @@ FE_ABF::interpolate( // TODO: Check what the face_orientation, face_flip and face_rotation have to be in 3D unsigned int k = QProjector::DataSetDescriptor::face (face, false, false, false, n_face_points); for (unsigned int i=0; i::unit_normal_direction[face]][k + fp]; } } @@ -774,6 +728,8 @@ FE_ABF::memory_consumption () const } -template class FE_ABF; + +/*-------------- Explicit Instantiations -------------------------------*/ +#include "fe_abf.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_abf.inst.in b/deal.II/source/fe/fe_abf.inst.in new file mode 100644 index 0000000000..39c7426305 --- /dev/null +++ b/deal.II/source/fe/fe_abf.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_ABF; + } + diff --git a/deal.II/source/fe/fe_bdm.cc b/deal.II/source/fe/fe_bdm.cc index 724e3779fe..be985b7cf4 100644 --- a/deal.II/source/fe/fe_bdm.cc +++ b/deal.II/source/fe/fe_bdm.cc @@ -27,6 +27,9 @@ #include #include + +DEAL_II_NAMESPACE_OPEN + template FE_BDM::FE_BDM (const unsigned int deg) : @@ -206,21 +209,6 @@ FE_BDM::interpolate( -#if deal_II_dimension == 1 - -template <> -std::vector -FE_BDM<1>::get_dpo_vector (const unsigned int deg) -{ - std::vector dpo(2); - dpo[0] = 1; - dpo[1] = deg; - return dpo; -} - -#endif - - template std::vector FE_BDM::get_dpo_vector (const unsigned int deg) @@ -233,7 +221,7 @@ FE_BDM::get_dpo_vector (const unsigned int deg) // and then there are interior dofs unsigned int - interior_dofs = dim*deg*(deg-1)/2; + interior_dofs = dim==1?deg:dim*deg*(deg-1)/2; if (dim>2) { interior_dofs *= deg-2; @@ -249,23 +237,16 @@ FE_BDM::get_dpo_vector (const unsigned int deg) -#if deal_II_dimension == 1 - -template <> -std::vector -FE_BDM<1>::get_ria_vector (const unsigned int) -{ - Assert (false, ExcImpossibleInDim(1)); - return std::vector(); -} - -#endif - - template std::vector FE_BDM::get_ria_vector (const unsigned int deg) { + if (dim==1) + { + Assert (false, ExcImpossibleInDim(1)); + return std::vector(); + } + Assert(dim==2, ExcNotImplemented()); const unsigned int dofs_per_cell = PolynomialsBDM::compute_n_pols(deg); const unsigned int dofs_per_face = deg+1; @@ -374,5 +355,8 @@ FE_BDM::initialize_support_points (const unsigned int deg) } -template class FE_BDM; +/*-------------- Explicit Instantiations -------------------------------*/ +#include "fe_bdm.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_bdm.inst.in b/deal.II/source/fe/fe_bdm.inst.in new file mode 100644 index 0000000000..314d27c5e9 --- /dev/null +++ b/deal.II/source/fe/fe_bdm.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_BDM; + } + diff --git a/deal.II/source/fe/fe_data.all_dimensions.cc b/deal.II/source/fe/fe_data.cc similarity index 100% rename from deal.II/source/fe/fe_data.all_dimensions.cc rename to deal.II/source/fe/fe_data.cc diff --git a/deal.II/source/fe/fe_dg_vector.cc b/deal.II/source/fe/fe_dg_vector.cc index 9c2c7d8f81..3cc8624b4c 100644 --- a/deal.II/source/fe/fe_dg_vector.cc +++ b/deal.II/source/fe/fe_dg_vector.cc @@ -18,10 +18,7 @@ DEAL_II_NAMESPACE_OPEN -template class FE_DGVector, deal_II_dimension>; -template class FE_DGVector, deal_II_dimension>; -template class FE_DGVector, deal_II_dimension>; -template class FE_DGVector, deal_II_dimension>; +#include "fe_dg_vector.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_dg_vector.inst.in b/deal.II/source/fe/fe_dg_vector.inst.in new file mode 100644 index 0000000000..7fbb81ab44 --- /dev/null +++ b/deal.II/source/fe/fe_dg_vector.inst.in @@ -0,0 +1,22 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_DGVector, deal_II_dimension>; + template class FE_DGVector, deal_II_dimension>; + template class FE_DGVector, deal_II_dimension>; + template class FE_DGVector, deal_II_dimension>; + } + diff --git a/deal.II/source/fe/fe_dgp.cc b/deal.II/source/fe/fe_dgp.cc index 9470a0bb65..802c3cd3c6 100644 --- a/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/source/fe/fe_dgp.cc @@ -53,7 +53,7 @@ FE_DGP::get_name () const // this function returns, so they // have to be kept in synch - std::ostringstream namebuf; + std::ostringstream namebuf; namebuf << "FE_DGP<" << dim << ">(" << this->degree << ")"; return namebuf.str(); @@ -110,7 +110,7 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, (dynamic_cast(&x_source_fe) != 0), typename FE:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -141,7 +141,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe (dynamic_cast(&x_source_fe) != 0), typename FE:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -257,10 +257,8 @@ FE_DGP::memory_consumption () const -template class FE_DGP; +// explicit instantiations +#include "fe_dgp.inst" -#if deal_II_dimension != 3 -template class FE_DGP; -#endif DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_dgp.inst.in b/deal.II/source/fe/fe_dgp.inst.in new file mode 100644 index 0000000000..053e3a3ca1 --- /dev/null +++ b/deal.II/source/fe/fe_dgp.inst.in @@ -0,0 +1,23 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_DGP; + +#if deal_II_dimension != 3 + template class FE_DGP; +#endif + } + diff --git a/deal.II/source/fe/fe_dgp_monomial.cc b/deal.II/source/fe/fe_dgp_monomial.cc index 03c88b2729..d7af3c4a8d 100644 --- a/deal.II/source/fe/fe_dgp_monomial.cc +++ b/deal.II/source/fe/fe_dgp_monomial.cc @@ -19,19 +19,18 @@ DEAL_II_NAMESPACE_OPEN // namespace for some functions that are used in this file. -namespace +namespace { // storage of hand-chosen support // points // -#if deal_II_dimension == 2 // For dim=2, dofs_per_cell of // FE_DGPMonomial(k) is given by // 0.5(k+1)(k+2), i.e. - // + // // k 0 1 2 3 4 5 6 7 // dofs 1 3 6 10 15 21 28 36 - // + // // indirect access of unit points: // the points for degree k are // located at @@ -45,13 +44,11 @@ namespace {0,0},{1,0},{0,1},{1,1},{1./3.,0},{2./3.,0},{0,1./3.},{0,2./3.},{0.5,1},{1,0.5}, {0,0},{1,0},{0,1},{1,1},{0.25,0},{0.5,0},{0.75,0},{0,0.25},{0,0.5},{0,0.75},{1./3.,1},{2./3.,1},{1,1./3.},{1,2./3.},{0.5,0.5} }; -#endif -#if deal_II_dimension == 3 // For dim=3, dofs_per_cell of // FE_DGPMonomial(k) is given by // 1./6.(k+1)(k+2)(k+3), i.e. - // + // // k 0 1 2 3 4 5 6 7 // dofs 1 4 10 20 35 56 84 120 const unsigned int start_index3d[6]={0,1,5,15/*,35*/}; @@ -60,14 +57,12 @@ namespace {0,0,0},{1,0,0},{0,1,0},{0,0,1}, {0,0,0},{1,0,0},{0,1,0},{0,0,1},{0.5,0,0},{0,0.5,0},{0,0,0.5},{1,1,0},{1,0,1},{0,1,1} }; -#endif - + template void generate_unit_points (const unsigned int, std::vector > &); -#if deal_II_dimension == 1 template <> void generate_unit_points (const unsigned int k, std::vector > &p) @@ -77,9 +72,7 @@ namespace for (unsigned int i=0; i void generate_unit_points (const unsigned int k, std::vector > &p) @@ -92,9 +85,7 @@ namespace p[i](1)=points2d[start_index2d[k]+i][1]; } } -#endif - -#if deal_II_dimension == 3 + template <> void generate_unit_points (const unsigned int k, std::vector > &p) @@ -108,7 +99,6 @@ namespace p[i](2)=points3d[start_index3d[k]+i][2]; } } -#endif } @@ -125,7 +115,7 @@ FE_DGPMonomial::FE_DGPMonomial (const unsigned int degree) { Assert(this->poly_space.n()==this->dofs_per_cell, ExcInternalError()); Assert(this->poly_space.degree()==this->degree, ExcInternalError()); - + // DG doesn't have constraints, so // leave them empty @@ -175,7 +165,7 @@ void FE_DGPMonomial:: get_interpolation_matrix (const FiniteElement &source_fe, FullMatrix &interpolation_matrix) const -{ +{ const FE_DGPMonomial *source_dgp_monomial = dynamic_cast *>(&source_fe); @@ -189,11 +179,11 @@ get_interpolation_matrix (const FiniteElement &source_fe, Assert (m == this->dofs_per_cell, ExcDimensionMismatch (m, this->dofs_per_cell)); Assert (n == source_dgp_monomial->dofs_per_cell, ExcDimensionMismatch (n, source_dgp_monomial->dofs_per_cell)); - + const unsigned int min_mn= interpolation_matrix.m() &source_fe, this_matrix(k,j)=this->poly_space.compute_value (j, unit_points[k]); this_matrix.gauss_jordan(); - + this_matrix.mmult(interpolation_matrix, source_fe_matrix); } } @@ -265,7 +255,7 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, (dynamic_cast*>(&x_source_fe) != 0), typename FiniteElement:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -294,7 +284,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, (dynamic_cast*>(&x_source_fe) != 0), typename FiniteElement:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -389,8 +379,6 @@ compare_for_face_domination (const FiniteElement &fe_other) const -#if deal_II_dimension == 1 - template <> bool FE_DGPMonomial<1>::has_support_on_face (const unsigned int, @@ -399,9 +387,7 @@ FE_DGPMonomial<1>::has_support_on_face (const unsigned int, return face_index==1 || (face_index==0 && this->degree==0); } -#endif -#if deal_II_dimension == 2 template <> bool @@ -422,9 +408,7 @@ FE_DGPMonomial<2>::has_support_on_face (const unsigned int shape_index, return support_on_face; } -#endif -#if deal_II_dimension == 3 template <> bool @@ -446,7 +430,6 @@ FE_DGPMonomial<3>::has_support_on_face (const unsigned int shape_index, return support_on_face; } -#endif template @@ -459,6 +442,8 @@ FE_DGPMonomial::memory_consumption () const -template class FE_DGPMonomial; +// explicit instantiations +#include "fe_dgp_monomial.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_dgp_monomial.inst.in b/deal.II/source/fe/fe_dgp_monomial.inst.in new file mode 100644 index 0000000000..401fe4df15 --- /dev/null +++ b/deal.II/source/fe/fe_dgp_monomial.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_DGPMonomial; + } + diff --git a/deal.II/source/fe/fe_dgp_nonparametric.cc b/deal.II/source/fe/fe_dgp_nonparametric.cc index 369d8ba1c6..1c2a80173a 100644 --- a/deal.II/source/fe/fe_dgp_nonparametric.cc +++ b/deal.II/source/fe/fe_dgp_nonparametric.cc @@ -24,7 +24,7 @@ DEAL_II_NAMESPACE_OPEN -template +template FE_DGPNonparametric::FE_DGPNonparametric (const unsigned int degree) : FiniteElement ( @@ -47,7 +47,7 @@ FE_DGPNonparametric::FE_DGPNonparametric (const unsigned int degre // refinement is not // implemented so far continue; - + const unsigned int nc = GeometryInfo::n_children(RefinementCase(ref_case)); for (unsigned int i=0;i::FE_DGPNonparametric (const unsigned int degre -template +template std::string FE_DGPNonparametric::get_name () const { @@ -109,7 +109,7 @@ FE_DGPNonparametric::get_name () const -template +template FiniteElement * FE_DGPNonparametric::clone() const { @@ -118,7 +118,7 @@ FE_DGPNonparametric::clone() const -template +template double FE_DGPNonparametric::shape_value (const unsigned int i, const Point &p) const @@ -129,7 +129,7 @@ FE_DGPNonparametric::shape_value (const unsigned int i, -template +template double FE_DGPNonparametric::shape_value_component (const unsigned int i, const Point &p, @@ -142,7 +142,7 @@ FE_DGPNonparametric::shape_value_component (const unsigned int i, -template +template Tensor<1,dim> FE_DGPNonparametric::shape_grad (const unsigned int i, const Point &p) const @@ -152,7 +152,7 @@ FE_DGPNonparametric::shape_grad (const unsigned int i, } -template +template Tensor<1,dim> FE_DGPNonparametric::shape_grad_component (const unsigned int i, const Point &p, @@ -165,7 +165,7 @@ FE_DGPNonparametric::shape_grad_component (const unsigned int i, -template +template Tensor<2,dim> FE_DGPNonparametric::shape_grad_grad (const unsigned int i, const Point &p) const @@ -176,7 +176,7 @@ FE_DGPNonparametric::shape_grad_grad (const unsigned int i, -template +template Tensor<2,dim> FE_DGPNonparametric::shape_grad_grad_component (const unsigned int i, const Point &p, @@ -193,7 +193,7 @@ FE_DGPNonparametric::shape_grad_grad_component (const unsigned int //--------------------------------------------------------------------------- -template +template std::vector FE_DGPNonparametric::get_dpo_vector (const unsigned int deg) { @@ -208,7 +208,7 @@ FE_DGPNonparametric::get_dpo_vector (const unsigned int deg) } -template +template UpdateFlags FE_DGPNonparametric::update_once (const UpdateFlags) const { @@ -221,7 +221,7 @@ FE_DGPNonparametric::update_once (const UpdateFlags) const } -template +template UpdateFlags FE_DGPNonparametric::update_each (const UpdateFlags flags) const { @@ -238,7 +238,7 @@ FE_DGPNonparametric::update_each (const UpdateFlags flags) const // Data field initialization //--------------------------------------------------------------------------- -template +template typename Mapping::InternalDataBase * FE_DGPNonparametric::get_data ( const UpdateFlags update_flags, @@ -283,7 +283,7 @@ FE_DGPNonparametric::get_data ( // Fill data of FEValues //--------------------------------------------------------------------------- -template +template void FE_DGPNonparametric::fill_fe_values ( const Mapping&, @@ -326,7 +326,7 @@ FE_DGPNonparametric::fill_fe_values ( -template +template void FE_DGPNonparametric::fill_fe_face_values ( const Mapping&, @@ -349,7 +349,7 @@ FE_DGPNonparametric::fill_fe_face_values ( Assert (flags & update_quadrature_points, ExcInternalError()); const unsigned int n_q_points = data.quadrature_points.size(); - + if (flags & (update_values | update_gradients)) for (unsigned int i=0; i::fill_fe_face_values ( -template +template void FE_DGPNonparametric::fill_fe_subface_values ( const Mapping&, @@ -388,12 +388,12 @@ FE_DGPNonparametric::fill_fe_subface_values ( Assert (dynamic_cast (&fedata) != 0, ExcInternalError()); InternalData &fe_data = static_cast (fedata); - + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); Assert (flags & update_quadrature_points, ExcInternalError()); const unsigned int n_q_points = data.quadrature_points.size(); - + if (flags & (update_values | update_gradients)) for (unsigned int i=0; i::fill_fe_subface_values ( -template +template unsigned int FE_DGPNonparametric::n_base_elements () const { @@ -422,7 +422,7 @@ FE_DGPNonparametric::n_base_elements () const -template +template const FiniteElement & FE_DGPNonparametric::base_element (const unsigned int index) const { @@ -432,7 +432,7 @@ FE_DGPNonparametric::base_element (const unsigned int index) const -template +template unsigned int FE_DGPNonparametric::element_multiplicity (const unsigned int index) const { @@ -442,7 +442,7 @@ FE_DGPNonparametric::element_multiplicity (const unsigned int inde -template +template void FE_DGPNonparametric:: get_face_interpolation_matrix (const FiniteElement &x_source_fe, @@ -460,7 +460,7 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, (dynamic_cast*>(&x_source_fe) != 0), typename FEE:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -471,7 +471,7 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, -template +template void FE_DGPNonparametric:: get_subface_interpolation_matrix (const FiniteElement &x_source_fe, @@ -490,7 +490,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe (dynamic_cast*>(&x_source_fe) != 0), typename FEE:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -501,7 +501,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe -template +template bool FE_DGPNonparametric::hp_constraints_are_implemented () const { @@ -510,7 +510,7 @@ FE_DGPNonparametric::hp_constraints_are_implemented () const -template +template std::vector > FE_DGPNonparametric:: hp_vertex_dof_identities (const FiniteElement &fe_other) const @@ -529,7 +529,7 @@ hp_vertex_dof_identities (const FiniteElement &fe_other) const -template +template std::vector > FE_DGPNonparametric:: hp_line_dof_identities (const FiniteElement &fe_other) const @@ -548,7 +548,7 @@ hp_line_dof_identities (const FiniteElement &fe_other) const -template +template std::vector > FE_DGPNonparametric:: hp_quad_dof_identities (const FiniteElement &fe_other) const @@ -567,7 +567,7 @@ hp_quad_dof_identities (const FiniteElement &fe_other) cons -template +template FiniteElementDomination::Domination FE_DGPNonparametric:: compare_for_face_domination (const FiniteElement &fe_other) const @@ -585,7 +585,7 @@ compare_for_face_domination (const FiniteElement &fe_other) const -template +template bool FE_DGPNonparametric::has_support_on_face (const unsigned int, const unsigned int) const @@ -595,7 +595,7 @@ FE_DGPNonparametric::has_support_on_face (const unsigned int, -template +template unsigned int FE_DGPNonparametric::memory_consumption () const { @@ -605,7 +605,7 @@ FE_DGPNonparametric::memory_consumption () const -template +template unsigned int FE_DGPNonparametric::get_degree () const { @@ -614,7 +614,8 @@ FE_DGPNonparametric::get_degree () const +// explicit instantiations +#include "fe_dgp_nonparametric.inst" -template class FE_DGPNonparametric; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_dgp_nonparametric.inst.in b/deal.II/source/fe/fe_dgp_nonparametric.inst.in new file mode 100644 index 0000000000..7e6cfe9486 --- /dev/null +++ b/deal.II/source/fe/fe_dgp_nonparametric.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_DGPNonparametric; + } + diff --git a/deal.II/source/fe/fe_dgq.cc b/deal.II/source/fe/fe_dgq.cc index 6afd3a0832..4a5c3192ef 100644 --- a/deal.II/source/fe/fe_dgq.cc +++ b/deal.II/source/fe/fe_dgq.cc @@ -752,14 +752,8 @@ FE_DGQArbitraryNodes::clone() const -template class FE_DGQ; +// explicit instantiations +#include "fe_dgq.inst" -//codimension 1 - -#if deal_II_dimension != 3 -template class FE_DGQ; -#endif - -template class FE_DGQArbitraryNodes; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_dgq.inst.in b/deal.II/source/fe/fe_dgq.inst.in new file mode 100644 index 0000000000..00c3243fae --- /dev/null +++ b/deal.II/source/fe/fe_dgq.inst.in @@ -0,0 +1,25 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_DGQ; + +#if deal_II_dimension != 3 + template class FE_DGQ; +#endif + + template class FE_DGQArbitraryNodes; + } + diff --git a/deal.II/source/fe/fe_face.cc b/deal.II/source/fe/fe_face.cc index 2d293c27ed..011c9dfe3a 100644 --- a/deal.II/source/fe/fe_face.cc +++ b/deal.II/source/fe/fe_face.cc @@ -46,7 +46,7 @@ FE_FaceQ::get_name () const // this function returns, so they // have to be kept in synch - std::ostringstream namebuf; + std::ostringstream namebuf; namebuf << "FE_FaceQ<" << dim << ">(" << this->degree << ")"; return namebuf.str(); @@ -62,8 +62,8 @@ FE_FaceQ::has_support_on_face ( { return (face_index == (shape_index/this->dofs_per_face)); } - - + + template std::vector FE_FaceQ::get_dpo_vector (const unsigned int deg) @@ -76,13 +76,9 @@ FE_FaceQ::get_dpo_vector (const unsigned int deg) } - -#if deal_II_dimension > 1 - template class FE_PolyFace >; -//template class FE_PolyFace, deal_II_dimension>; -//template class FE_PolyFace, deal_II_dimension>; - template class FE_FaceQ; -#endif - +// explicit instantiations +#include "fe_face.inst" + + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_face.inst.in b/deal.II/source/fe/fe_face.inst.in new file mode 100644 index 0000000000..631a8f52f6 --- /dev/null +++ b/deal.II/source/fe/fe_face.inst.in @@ -0,0 +1,24 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { +#if deal_II_dimension > 1 + template class FE_PolyFace >; + //template class FE_PolyFace, deal_II_dimension>; + //template class FE_PolyFace, deal_II_dimension>; + template class FE_FaceQ; +#endif + } + diff --git a/deal.II/source/fe/fe_nedelec.cc b/deal.II/source/fe/fe_nedelec.cc index 47b4c721f5..e904eb39c0 100644 --- a/deal.II/source/fe/fe_nedelec.cc +++ b/deal.II/source/fe/fe_nedelec.cc @@ -173,7 +173,6 @@ FiniteElement //--------------------------------------------------------------------------- -#if deal_II_dimension == 1 // Set the generalized support // points and precompute the @@ -188,14 +187,14 @@ FE_Nedelec<1>::initialize_support_points (const unsigned int) Assert (false, ExcNotImplemented ()); } -#elif deal_II_dimension == 2 + template <> void FE_Nedelec<2>::initialize_support_points (const unsigned int degree) { const int dim = 2; - + // Create polynomial basis. const std::vector >& lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis (degree + 1); @@ -280,14 +279,14 @@ FE_Nedelec<2>::initialize_support_points (const unsigned int degree) } } -#else + template <> void FE_Nedelec<3>::initialize_support_points (const unsigned int degree) { const int dim = 3; - + // Create polynomial basis. const std::vector >& lobatto_polynomials = Polynomials::Lobatto::generate_complete_basis (degree + 1); @@ -460,11 +459,8 @@ FE_Nedelec<3>::initialize_support_points (const unsigned int degree) } } -#endif -#if deal_II_dimension == 1 - // Set the restriction matrices. template <> void @@ -476,7 +472,7 @@ FE_Nedelec<1>::initialize_restriction () this->restriction[0][i].reinit(0, 0); } -#endif + // Restriction operator template @@ -496,7 +492,7 @@ FE_Nedelec::initialize_restriction () n_edge_quadrature_points = edge_quadrature.size (); const unsigned int index = RefinementCase::isotropic_refinement - 1; - + switch (dim) { case 2: @@ -510,12 +506,12 @@ FE_Nedelec::initialize_restriction () ++q_point) { const double weight = 2.0 * edge_quadrature.weight (q_point); - + if (edge_quadrature_points[q_point] (0) < 0.5) { Point quadrature_point (0.0, 2.0 * edge_quadrature_points[q_point] (0)); - + this->restriction[index][0] (0, dof) += weight * this->shape_value_component (dof, @@ -538,13 +534,13 @@ FE_Nedelec::initialize_restriction () quadrature_point, 0); } - + else { Point quadrature_point (0.0, 2.0 * edge_quadrature_points[q_point] (0) - 1.0); - + this->restriction[index][2] (0, dof) += weight * this->shape_value_component (dof, @@ -568,7 +564,7 @@ FE_Nedelec::initialize_restriction () 0); } } - + // Then project the shape functions // of the child cells to the higher // order shape functions of the @@ -579,38 +575,38 @@ FE_Nedelec::initialize_restriction () legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg); FullMatrix system_matrix_inv (deg, deg); - + { FullMatrix assembling_matrix (deg, n_edge_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = std::sqrt (edge_quadrature.weight (q_point)); - + for (unsigned int i = 0; i < deg; ++i) assembling_matrix (i, q_point) = weight * legendre_polynomials[i + 1].value (edge_quadrature_points[q_point] (0)); } - + FullMatrix system_matrix (deg, deg); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.invert (system_matrix); } - + FullMatrix solution (deg, 4); FullMatrix system_rhs (deg, 4); Vector tmp (4); - + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) for (unsigned int i = 0; i < 2; ++i) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { @@ -621,12 +617,12 @@ FE_Nedelec::initialize_restriction () const Point quadrature_point_1 (edge_quadrature_points[q_point] (0), i); - + if (edge_quadrature_points[q_point] (0) < 0.5) { Point quadrature_point_2 (i, 2.0 * edge_quadrature_points[q_point] (0)); - + tmp (0) = weight * (2.0 * this->shape_value_component (dof, quadrature_point_2, 1) @@ -659,7 +655,7 @@ FE_Nedelec::initialize_restriction () ((i + 2) * this->degree, quadrature_point_1, 0); } - + else { tmp (0) = -1.0 * weight @@ -668,11 +664,11 @@ FE_Nedelec::initialize_restriction () * this->shape_value_component (i * this->degree, quadrature_point_0, 1); - + Point quadrature_point_2 (i, 2.0 * edge_quadrature_points[q_point] (0) - 1.0); - + tmp (1) = weight * (2.0 * this->shape_value_component (dof, quadrature_point_2, 1) @@ -699,20 +695,20 @@ FE_Nedelec::initialize_restriction () ((i + 2) * this->degree, quadrature_point_1, 0)); } - + for (unsigned int j = 0; j < deg; ++j) { const double L_j = legendre_polynomials[j + 1].value (edge_quadrature_points[q_point] (0)); - + for (unsigned int k = 0; k < tmp.size (); ++k) system_rhs (j, k) += tmp (k) * L_j; } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int j = 0; j < deg; ++j) for (unsigned int k = 0; k < 2; ++k) { @@ -720,14 +716,14 @@ FE_Nedelec::initialize_restriction () this->restriction[index][i + 2 * k] (i * this->degree + j + 1, dof) = solution (j, k); - + if (std::abs (solution (j, k + 2)) > 1e-14) this->restriction[index][2 * i + k] ((i + 2) * this->degree + j + 1, dof) = solution (j, k + 2); } } - + const QGauss quadrature (2 * this->degree); const std::vector >& quadrature_points = quadrature.get_points (); @@ -738,51 +734,51 @@ FE_Nedelec::initialize_restriction () const unsigned int n_boundary_dofs = GeometryInfo::faces_per_cell * this->degree; const unsigned int& n_quadrature_points = quadrature.size (); - + { FullMatrix assembling_matrix (deg * this->degree, n_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { const double weight = std::sqrt (quadrature.weight (q_point)); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0)); - + for (unsigned int j = 0; j < deg; ++j) assembling_matrix (i * deg + j, q_point) = L_i * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1)); } } - + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ()); system_matrix_inv.invert (system_matrix); } - + solution.reinit (system_matrix_inv.m (), 8); system_rhs.reinit (system_matrix_inv.m (), 8); tmp.reinit (8); - + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { tmp = 0.0; - + if (quadrature_points[q_point] (0) < 0.5) { if (quadrature_points[q_point] (1) < 0.5) @@ -790,27 +786,27 @@ FE_Nedelec::initialize_restriction () const Point quadrature_point (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1)); - + tmp (0) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (1) += 2.0 * this->shape_value_component (dof, quadrature_point, 1); } - + else { const Point quadrature_point (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1) - 1.0); - + tmp (4) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (5) += 2.0 * this->shape_value_component (dof, quadrature_point, 1); } } - + else if (quadrature_points[q_point] (1) < 0.5) { @@ -818,13 +814,13 @@ FE_Nedelec::initialize_restriction () (2.0 * quadrature_points[q_point] (0) - 1.0, 2.0 * quadrature_points[q_point] (1)); - + tmp (2) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (3) += 2.0 * this->shape_value_component (dof, quadrature_point, 1); } - + else { const Point quadrature_point @@ -832,13 +828,13 @@ FE_Nedelec::initialize_restriction () - 1.0, 2.0 * quadrature_points[q_point] (1) - 1.0); - + tmp (6) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (7) += 2.0 * this->shape_value_component (dof, quadrature_point, 1); } - + for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j <= deg; ++j) { @@ -864,9 +860,9 @@ FE_Nedelec::initialize_restriction () (i * this->degree + j, quadrature_points[q_point], 1); } - + tmp *= quadrature.weight (q_point); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i_0 @@ -875,7 +871,7 @@ FE_Nedelec::initialize_restriction () const double L_i_1 = legendre_polynomials[i].value (quadrature_points[q_point] (1)); - + for (unsigned int j = 0; j < deg; ++j) { const double l_j_0 @@ -884,7 +880,7 @@ FE_Nedelec::initialize_restriction () const double l_j_1 = L_i_1 * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0)); - + for (unsigned int k = 0; k < 4; ++k) { system_rhs (i * deg + j, 2 * k) @@ -895,9 +891,9 @@ FE_Nedelec::initialize_restriction () } } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) for (unsigned int k = 0; k < 4; ++k) @@ -907,7 +903,7 @@ FE_Nedelec::initialize_restriction () this->restriction[index][k] (i * deg + j + n_boundary_dofs, dof) = solution (i * deg + j, 2 * k); - + if (std::abs (solution (i * deg + j, 2 * k + 1)) > 1e-14) this->restriction[index][k] @@ -932,7 +928,7 @@ FE_Nedelec::initialize_restriction () ++q_point) { const double weight = 2.0 * edge_quadrature.weight (q_point); - + if (edge_quadrature_points[q_point] (0) < 0.5) for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) @@ -940,7 +936,7 @@ FE_Nedelec::initialize_restriction () Point quadrature_point (i, 2.0 * edge_quadrature_points[q_point] (0), j); - + this->restriction[index][i + 4 * j] ((i + 4 * j) * this->degree, dof) += weight * this->shape_value_component (dof, @@ -962,7 +958,7 @@ FE_Nedelec::initialize_restriction () quadrature_point, 2); } - + else for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) @@ -970,7 +966,7 @@ FE_Nedelec::initialize_restriction () Point quadrature_point (i, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, j); - + this->restriction[index][i + 4 * j + 2] ((i + 4 * j) * this->degree, dof) += weight * this->shape_value_component (dof, @@ -994,7 +990,7 @@ FE_Nedelec::initialize_restriction () 2); } } - + // Then project the shape functions // of the child cells to the higher // order shape functions of the @@ -1005,39 +1001,39 @@ FE_Nedelec::initialize_restriction () legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg); FullMatrix system_matrix_inv (deg, deg); - + { FullMatrix assembling_matrix (deg, n_edge_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = std::sqrt (edge_quadrature.weight (q_point)); - + for (unsigned int i = 0; i < deg; ++i) assembling_matrix (i, q_point) = weight * legendre_polynomials[i + 1].value (edge_quadrature_points[q_point] (0)); } - + FullMatrix system_matrix (deg, deg); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.invert (system_matrix); } - + FullMatrix solution (deg, 6); FullMatrix system_rhs (deg, 6); Vector tmp (6); - + for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { @@ -1051,13 +1047,13 @@ FE_Nedelec::initialize_restriction () (edge_quadrature_points[q_point] (0), i, j); const Point quadrature_point_2 (i, j, edge_quadrature_points[q_point] (0)); - + if (edge_quadrature_points[q_point] (0) < 0.5) { Point quadrature_point_3 (i, 2.0 * edge_quadrature_points[q_point] (0), j); - + tmp (0) = weight * (2.0 * this->shape_value_component (dof, quadrature_point_3, 1) @@ -1112,7 +1108,7 @@ FE_Nedelec::initialize_restriction () ((i + 2 * (j + 4)) * this->degree, quadrature_point_2, 2); } - + else { tmp (0) = -1.0 * weight @@ -1122,11 +1118,11 @@ FE_Nedelec::initialize_restriction () * this->shape_value_component ((i + 4 * j) * this->degree, quadrature_point_0, 1); - + Point quadrature_point_3 (i, 2.0 * edge_quadrature_points[q_point] (0) - 1.0, j); - + tmp (1) = weight * (2.0 * this->shape_value_component (dof, quadrature_point_3, 1) @@ -1175,20 +1171,20 @@ FE_Nedelec::initialize_restriction () ((i + 2 * (j + 4)) * this->degree, quadrature_point_2, 2)); } - + for (unsigned int k = 0; k < deg; ++k) { const double L_k = legendre_polynomials[k + 1].value (edge_quadrature_points[q_point] (0)); - + for (unsigned int l = 0; l < tmp.size (); ++l) system_rhs (k, l) += tmp (l) * L_k; } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int k = 0; k < 2; ++k) for (unsigned int l = 0; l < deg; ++l) { @@ -1196,12 +1192,12 @@ FE_Nedelec::initialize_restriction () this->restriction[index][i + 2 * (2 * j + k)] ((i + 4 * j) * this->degree + l + 1, dof) = solution (l, k); - + if (std::abs (solution (l, k + 2)) > 1e-14) this->restriction[index][2 * (i + 2 * j) + k] ((i + 4 * j + 2) * this->degree + l + 1, dof) = solution (l, k + 2); - + if (std::abs (solution (l, k + 4)) > 1e-14) this->restriction[index][i + 2 * (j + 2 * k)] ((i + 2 * (j + 4)) * this->degree + l + 1, @@ -1209,7 +1205,7 @@ FE_Nedelec::initialize_restriction () = solution (l, k + 4); } } - + const QGauss<2> face_quadrature (2 * this->degree); const std::vector >& face_quadrature_points = face_quadrature.get_points (); @@ -1221,55 +1217,55 @@ FE_Nedelec::initialize_restriction () = GeometryInfo::lines_per_cell * this->degree; const unsigned int& n_face_quadrature_points = face_quadrature.size (); - + { FullMatrix assembling_matrix (deg * this->degree, n_face_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_face_quadrature_points; ++q_point) { const double weight = std::sqrt (face_quadrature.weight (q_point)); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i = weight * legendre_polynomials[i].value (face_quadrature_points[q_point] (0)); - + for (unsigned int j = 0; j < deg; ++j) assembling_matrix (i * deg + j, q_point) = L_i * lobatto_polynomials[j + 2].value (face_quadrature_points[q_point] (1)); } } - + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ()); system_matrix_inv.invert (system_matrix); } - + solution.reinit (system_matrix_inv.m (), 24); system_rhs.reinit (system_matrix_inv.m (), 24); tmp.reinit (24); - + for (unsigned int i = 0; i < 2; ++i) for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_face_quadrature_points; ++q_point) { tmp = 0.0; - + if (face_quadrature_points[q_point] (0) < 0.5) { if (face_quadrature_points[q_point] (1) < 0.5) @@ -1277,7 +1273,7 @@ FE_Nedelec::initialize_restriction () Point quadrature_point_0 (i, 2.0 * face_quadrature_points[q_point] (0), 2.0 * face_quadrature_points[q_point] (1)); - + tmp (0) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); tmp (1) += 2.0 * this->shape_value_component @@ -1299,14 +1295,14 @@ FE_Nedelec::initialize_restriction () tmp (17) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); } - + else { Point quadrature_point_0 (i, 2.0 * face_quadrature_points[q_point] (0), 2.0 * face_quadrature_points[q_point] (1) - 1.0); - + tmp (2) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); tmp (3) += 2.0 * this->shape_value_component @@ -1330,7 +1326,7 @@ FE_Nedelec::initialize_restriction () (dof, quadrature_point_0, 1); } } - + else if (face_quadrature_points[q_point] (1) < 0.5) { @@ -1338,7 +1334,7 @@ FE_Nedelec::initialize_restriction () 2.0 * face_quadrature_points[q_point] (0) - 1.0, 2.0 * face_quadrature_points[q_point] (1)); - + tmp (4) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); tmp (5) += 2.0 * this->shape_value_component @@ -1361,7 +1357,7 @@ FE_Nedelec::initialize_restriction () tmp (21) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); } - + else { Point quadrature_point_0 (i, @@ -1369,7 +1365,7 @@ FE_Nedelec::initialize_restriction () - 1.0, 2.0 * face_quadrature_points[q_point] (1) - 1.0); - + tmp (6) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); tmp (7) += 2.0 * this->shape_value_component @@ -1393,7 +1389,7 @@ FE_Nedelec::initialize_restriction () tmp (23) += 2.0 * this->shape_value_component (dof, quadrature_point_0, 1); } - + const Point quadrature_point_0 (i, face_quadrature_points[q_point] (0), face_quadrature_points[q_point] (1)); @@ -1405,7 +1401,7 @@ FE_Nedelec::initialize_restriction () (face_quadrature_points[q_point] (0), face_quadrature_points[q_point] (1), i); - + for (unsigned int j = 0; j < 2; ++j) for (unsigned int k = 0; k < 2; ++k) for (unsigned int l = 0; l <= deg; ++l) @@ -1451,9 +1447,9 @@ FE_Nedelec::initialize_restriction () ((4 * i + k) * this->degree + l, quadrature_point_2, 1); } - + tmp *= face_quadrature.weight (q_point); - + for (unsigned int j = 0; j <= deg; ++j) { const double L_j_0 @@ -1462,7 +1458,7 @@ FE_Nedelec::initialize_restriction () const double L_j_1 = legendre_polynomials[j].value (face_quadrature_points[q_point] (1)); - + for (unsigned int k = 0; k < deg; ++k) { const double l_k_0 @@ -1471,7 +1467,7 @@ FE_Nedelec::initialize_restriction () const double l_k_1 = L_j_1 * lobatto_polynomials[k + 2].value (face_quadrature_points[q_point] (0)); - + for (unsigned int l = 0; l < 4; ++l) { system_rhs (j * deg + k, 2 * l) @@ -1490,9 +1486,9 @@ FE_Nedelec::initialize_restriction () } } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int j = 0; j < 2; ++j) for (unsigned int k = 0; k < 2; ++k) for (unsigned int l = 0; l <= deg; ++l) @@ -1503,10 +1499,10 @@ FE_Nedelec::initialize_restriction () > 1e-14) this->restriction[index][i + 2 * (2 * j + k)] ((2 * i * this->degree + l) * deg + m - + n_edge_dofs, + + n_edge_dofs, dof) = solution (l * deg + m, 2 * (j + 2 * k)); - + if (std::abs (solution (l * deg + m, 2 * (j + 2 * k) + 1)) > 1e-14) @@ -1515,7 +1511,7 @@ FE_Nedelec::initialize_restriction () + n_edge_dofs, dof) = solution (l * deg + m, 2 * (j + 2 * k) + 1); - + if (std::abs (solution (l * deg + m, 2 * (j + 2 * (k + 2)))) > 1e-14) @@ -1524,7 +1520,7 @@ FE_Nedelec::initialize_restriction () + n_edge_dofs, dof) = solution (l * deg + m, 2 * (j + 2 * (k + 2))); - + if (std::abs (solution (l * deg + m, 2 * (j + 2 * k) + 9)) > 1e-14) @@ -1533,7 +1529,7 @@ FE_Nedelec::initialize_restriction () + n_edge_dofs, dof) = solution (l * deg + m, 2 * (j + 2 * k) + 9); - + if (std::abs (solution (l * deg + m, 2 * (j + 2 * (k + 4)))) > 1e-14) @@ -1542,7 +1538,7 @@ FE_Nedelec::initialize_restriction () + n_edge_dofs, dof) = solution (l * deg + m, 2 * (j + 2 * (k + 4))); - + if (std::abs (solution (l * deg + m, 2 * (j + 2 * k) + 17)) > 1e-14) @@ -1553,7 +1549,7 @@ FE_Nedelec::initialize_restriction () 2 * (j + 2 * k) + 17); } } - + const QGauss quadrature (2 * this->degree); const std::vector >& quadrature_points = quadrature.get_points (); @@ -1561,30 +1557,30 @@ FE_Nedelec::initialize_restriction () = 2 * GeometryInfo::faces_per_cell * deg * this->degree + n_edge_dofs; const unsigned int& n_quadrature_points = quadrature.size (); - + { FullMatrix assembling_matrix (deg * deg * this->degree, n_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { const double weight = std::sqrt (quadrature.weight (q_point)); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0)); - + for (unsigned int j = 0; j < deg; ++j) { const double l_j = L_i * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1)); - + for (unsigned int k = 0; k < deg; ++k) assembling_matrix ((i * deg + j) * deg + k, q_point) @@ -1593,29 +1589,29 @@ FE_Nedelec::initialize_restriction () } } } - + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ()); system_matrix_inv.invert (system_matrix); } - + solution.reinit (system_matrix_inv.m (), 24); system_rhs.reinit (system_matrix_inv.m (), 24); tmp.reinit (24); - + for (unsigned int dof = 0; dof < this->dofs_per_cell; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { tmp = 0.0; - + if (quadrature_points[q_point] (0) < 0.5) { if (quadrature_points[q_point] (1) < 0.5) @@ -1626,7 +1622,7 @@ FE_Nedelec::initialize_restriction () (2.0 * quadrature_points[q_point] (0), 2.0 * quadrature_points[q_point] (1), 2.0 * quadrature_points[q_point] (2)); - + tmp (0) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (1) += 2.0 * this->shape_value_component @@ -1634,7 +1630,7 @@ FE_Nedelec::initialize_restriction () tmp (2) += 2.0 * this->shape_value_component (dof, quadrature_point, 2); } - + else { const Point quadrature_point @@ -1642,7 +1638,7 @@ FE_Nedelec::initialize_restriction () 2.0 * quadrature_points[q_point] (1), 2.0 * quadrature_points[q_point] (2) - 1.0); - + tmp (3) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (4) += 2.0 * this->shape_value_component @@ -1651,7 +1647,7 @@ FE_Nedelec::initialize_restriction () (dof, quadrature_point, 2); } } - + else if (quadrature_points[q_point] (2) < 0.5) { @@ -1660,7 +1656,7 @@ FE_Nedelec::initialize_restriction () 2.0 * quadrature_points[q_point] (1) - 1.0, 2.0 * quadrature_points[q_point] (2)); - + tmp (6) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (7) += 2.0 * this->shape_value_component @@ -1668,7 +1664,7 @@ FE_Nedelec::initialize_restriction () tmp (8) += 2.0 * this->shape_value_component (dof, quadrature_point, 2); } - + else { const Point quadrature_point @@ -1677,7 +1673,7 @@ FE_Nedelec::initialize_restriction () - 1.0, 2.0 * quadrature_points[q_point] (2) - 1.0); - + tmp (9) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (10) += 2.0 * this->shape_value_component @@ -1686,7 +1682,7 @@ FE_Nedelec::initialize_restriction () (dof, quadrature_point, 2); } } - + else if (quadrature_points[q_point] (1) < 0.5) { @@ -1697,7 +1693,7 @@ FE_Nedelec::initialize_restriction () - 1.0, 2.0 * quadrature_points[q_point] (1), 2.0 * quadrature_points[q_point] (2)); - + tmp (12) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (13) += 2.0 * this->shape_value_component @@ -1705,7 +1701,7 @@ FE_Nedelec::initialize_restriction () tmp (14) += 2.0 * this->shape_value_component (dof, quadrature_point, 2); } - + else { const Point quadrature_point @@ -1714,7 +1710,7 @@ FE_Nedelec::initialize_restriction () 2.0 * quadrature_points[q_point] (1), 2.0 * quadrature_points[q_point] (2) - 1.0); - + tmp (15) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (16) += 2.0 * this->shape_value_component @@ -1723,7 +1719,7 @@ FE_Nedelec::initialize_restriction () (dof, quadrature_point, 2); } } - + else if (quadrature_points[q_point] (2) < 0.5) { @@ -1733,7 +1729,7 @@ FE_Nedelec::initialize_restriction () 2.0 * quadrature_points[q_point] (1) - 1.0, 2.0 * quadrature_points[q_point] (2)); - + tmp (18) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (19) += 2.0 * this->shape_value_component @@ -1741,7 +1737,7 @@ FE_Nedelec::initialize_restriction () tmp (20) += 2.0 * this->shape_value_component (dof, quadrature_point, 2); } - + else { const Point quadrature_point @@ -1751,7 +1747,7 @@ FE_Nedelec::initialize_restriction () - 1.0, 2.0 * quadrature_points[q_point] (2) - 1.0); - + tmp (21) += 2.0 * this->shape_value_component (dof, quadrature_point, 0); tmp (22) += 2.0 * this->shape_value_component @@ -1759,7 +1755,7 @@ FE_Nedelec::initialize_restriction () tmp (23) += 2.0 * this->shape_value_component (dof, quadrature_point, 2); } - + for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) for (unsigned int k = 0; k < 2; ++k) @@ -1784,7 +1780,7 @@ FE_Nedelec::initialize_restriction () * this->shape_value_component ((2 * (j + 4) + k) * this->degree + l, quadrature_points[q_point], 2); - + for (unsigned int m = 0; m < deg; ++m) { tmp (3 * (i + 2 * (j + 2 * k))) @@ -1841,9 +1837,9 @@ FE_Nedelec::initialize_restriction () quadrature_points[q_point], 2); } } - + tmp *= quadrature.weight (q_point); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i_0 @@ -1855,7 +1851,7 @@ FE_Nedelec::initialize_restriction () const double L_i_2 = legendre_polynomials[i].value (quadrature_points[q_point] (2)); - + for (unsigned int j = 0; j < deg; ++j) { const double l_j_0 @@ -1867,7 +1863,7 @@ FE_Nedelec::initialize_restriction () const double l_j_2 = L_i_2 * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (0)); - + for (unsigned int k = 0; k < deg; ++k) { const double l_k_0 @@ -1879,7 +1875,7 @@ FE_Nedelec::initialize_restriction () const double l_k_2 = l_j_2 * lobatto_polynomials[k + 2].value (quadrature_points[q_point] (1)); - + for (unsigned int l = 0; l < 8; ++l) { system_rhs ((i * deg + j) * deg + k, @@ -1896,9 +1892,9 @@ FE_Nedelec::initialize_restriction () } } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j < 2; ++j) for (unsigned int k = 0; k < 2; ++k) @@ -1914,7 +1910,7 @@ FE_Nedelec::initialize_restriction () ((l * deg + m) * deg + n + n_boundary_dofs, dof) = solution ((l * deg + m) * deg + n, 3 * (i + 2 * (j + 2 * k))); - + if (std::abs (solution ((l * deg + m) * deg + n, 3 * (i + 2 * (j + 2 * k)) + 1)) @@ -1924,7 +1920,7 @@ FE_Nedelec::initialize_restriction () + n_boundary_dofs, dof) = solution ((l * deg + m) * deg + n, 3 * (i + 2 * (j + 2 * k)) + 1); - + if (std::abs (solution ((l * deg + m) * deg + n, 3 * (i + 2 * (j + 2 * k)) + 2)) @@ -1947,21 +1943,6 @@ FE_Nedelec::initialize_restriction () } -#if deal_II_dimension == 1 - -template <> -std::vector -FE_Nedelec<1>::get_dpo_vector (const unsigned int degree) -{ - std::vector dpo (2); - - dpo[0] = 1; - dpo[1] = degree; - return dpo; -} - -#endif - template std::vector @@ -2511,19 +2492,19 @@ const } } -#if deal_II_dimension == 1 -template + +template <> void -FE_Nedelec::get_subface_interpolation_matrix( - const FiniteElement&, +FE_Nedelec<1>::get_subface_interpolation_matrix( + const FiniteElement<1>&, const unsigned int, FullMatrix&) const { Assert (false, ExcNotImplemented ()); } -#else + // In this function we compute the // subface interpolation matrix. @@ -2589,7 +2570,7 @@ FE_Nedelec::get_subface_interpolation_matrix( const std::vector >& edge_quadrature_points = edge_quadrature.get_points (); const unsigned int& n_edge_quadrature_points = edge_quadrature.size (); - + switch (dim) { case 2: @@ -2601,49 +2582,49 @@ FE_Nedelec::get_subface_interpolation_matrix( const Point quadrature_point (0.0, 0.5 * (edge_quadrature_points[q_point] (0) + subface)); - + interpolation_matrix (0, dof) += 0.5 * edge_quadrature.weight (q_point) * this->shape_value_component (dof, quadrature_point, 1); } - + if (deg > 0) { const std::vector >& legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg); FullMatrix system_matrix_inv (deg, deg); - + { FullMatrix assembling_matrix (deg, n_edge_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = std::sqrt (edge_quadrature.weight (q_point)); - + for (unsigned int i = 0; i < deg; ++i) assembling_matrix (i, q_point) = weight * legendre_polynomials[i + 1].value (edge_quadrature_points[q_point] (0)); } - + FullMatrix system_matrix (deg, deg); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.invert (system_matrix); } - + Vector solution (deg); Vector system_rhs (deg); - + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { @@ -2659,21 +2640,21 @@ FE_Nedelec::get_subface_interpolation_matrix( dof) * source_fe.shape_value_component (0, quadrature_point_1, 1)); - + for (unsigned int i = 0; i < deg; ++i) system_rhs (i) += tmp * legendre_polynomials[i + 1].value (edge_quadrature_points[q_point] (0)); } - + system_matrix_inv.vmult (solution, system_rhs); - + for (unsigned int i = 0; i < deg; ++i) if (std::abs (solution (i)) > 1e-14) interpolation_matrix (i + 1, dof) = solution (i); } } - + break; } @@ -2681,20 +2662,20 @@ FE_Nedelec::get_subface_interpolation_matrix( { const double shifts[4][2] = { { 0.0, 0.0 }, { 1.0, 0.0 }, { 0.0, 1.0 }, { 1.0, 1.0 } }; - + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = 0.5 * edge_quadrature.weight (q_point); - + for (unsigned int i = 0; i < 2; ++i) { Point quadrature_point (0.5 * (i + shifts[subface][0]), 0.5 * (edge_quadrature_points[q_point] (0) + shifts[subface][1]), 0.0); - + interpolation_matrix (i * source_fe.degree, dof) += weight * this->shape_value_component (this->face_to_cell_index (dof, 4), @@ -2710,51 +2691,51 @@ FE_Nedelec::get_subface_interpolation_matrix( quadrature_point, 0); } } - + if (deg > 0) { const std::vector >& legendre_polynomials = Polynomials::Legendre::generate_complete_basis (deg); FullMatrix system_matrix_inv (deg, deg); - + { FullMatrix assembling_matrix (deg, n_edge_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = std::sqrt (edge_quadrature.weight (q_point)); - + for (unsigned int i = 0; i < deg; ++i) assembling_matrix (i, q_point) = weight * legendre_polynomials[i + 1].value (edge_quadrature_points[q_point] (0)); } - + FullMatrix system_matrix (deg, deg); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.invert (system_matrix); } - + FullMatrix solution (deg, GeometryInfo::lines_per_face); FullMatrix system_rhs (deg, GeometryInfo::lines_per_face); Vector tmp (GeometryInfo::lines_per_face); - + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_edge_quadrature_points; ++q_point) { const double weight = edge_quadrature.weight (q_point); - + for (unsigned int i = 0; i < 2; ++i) { Point @@ -2765,7 +2746,7 @@ FE_Nedelec::get_subface_interpolation_matrix( Point quadrature_point_1 (i, edge_quadrature_points[q_point] (0), 0.0); - + tmp (i) = weight * (0.5 * this->shape_value_component (this->face_to_cell_index (dof, 4), @@ -2794,21 +2775,21 @@ FE_Nedelec::get_subface_interpolation_matrix( ((i + 2) * source_fe.degree, quadrature_point_1, 0)); } - + for (unsigned int i = 0; i < deg; ++i) { const double L_i = legendre_polynomials[i + 1].value (edge_quadrature_points[q_point] (0)); - + for (unsigned int j = 0; j < GeometryInfo::lines_per_face; ++j) system_rhs (i, j) += tmp (j) * L_i; } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int i = 0; i < GeometryInfo::lines_per_face; ++i) for (unsigned int j = 0; j < deg; ++j) @@ -2816,7 +2797,7 @@ FE_Nedelec::get_subface_interpolation_matrix( interpolation_matrix (i * source_fe.degree + j + 1, dof) = solution (j, i); } - + const QGauss<2> quadrature (source_fe.degree); const std::vector >& quadrature_points = quadrature.get_points (); @@ -2827,46 +2808,46 @@ FE_Nedelec::get_subface_interpolation_matrix( const unsigned int n_boundary_dofs = GeometryInfo::lines_per_face * source_fe.degree; const unsigned int& n_quadrature_points = quadrature.size (); - + { FullMatrix assembling_matrix (deg * this->degree, n_quadrature_points); - + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { const double weight = quadrature.weight (q_point); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i = weight * legendre_polynomials[i].value (quadrature_points[q_point] (0)); - + for (unsigned int j = 0; j < deg; ++j) assembling_matrix (i * deg + j, q_point) = L_i * lobatto_polynomials[j + 2].value (quadrature_points[q_point] (1)); } } - + FullMatrix system_matrix (assembling_matrix.m (), assembling_matrix.m ()); - + assembling_matrix.mTmult (system_matrix, assembling_matrix); system_matrix_inv.reinit (system_matrix.m (), system_matrix.m ()); system_matrix_inv.invert (system_matrix); } - + solution.reinit (system_matrix_inv.m (), 2); system_rhs.reinit (system_matrix_inv.m (), 2); tmp.reinit (2); - + for (unsigned int dof = 0; dof < this->dofs_per_face; ++dof) { system_rhs = 0.0; - + for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { @@ -2885,7 +2866,7 @@ FE_Nedelec::get_subface_interpolation_matrix( quadrature_point = Point (quadrature_points[q_point] (0), quadrature_points[q_point] (1), 0.0); - + for (unsigned int i = 0; i < 2; ++i) for (unsigned int j = 0; j <= deg; ++j) { @@ -2900,16 +2881,16 @@ FE_Nedelec::get_subface_interpolation_matrix( (i * source_fe.degree + j, quadrature_point, 1); } - + tmp *= quadrature.weight (q_point); - + for (unsigned int i = 0; i <= deg; ++i) { const double L_i_0 = legendre_polynomials[i].value (quadrature_points[q_point] (0)); const double L_i_1 = legendre_polynomials[i].value (quadrature_points[q_point] (1)); - + for (unsigned int j = 0; j < deg; ++j) { system_rhs (i * deg + j, 0) += tmp (0) * L_i_0 @@ -2921,9 +2902,9 @@ FE_Nedelec::get_subface_interpolation_matrix( } } } - + system_matrix_inv.mmult (solution, system_rhs); - + for (unsigned int i = 0; i <= deg; ++i) for (unsigned int j = 0; j < deg; ++j) { @@ -2931,7 +2912,7 @@ FE_Nedelec::get_subface_interpolation_matrix( interpolation_matrix (i * deg + j + n_boundary_dofs, dof) = solution (i * deg + j, 0); - + if (std::abs (solution (i * deg + j, 1)) > 1e-14) interpolation_matrix (i + (j + deg) * source_fe.degree @@ -2940,7 +2921,7 @@ FE_Nedelec::get_subface_interpolation_matrix( } } } - + break; } @@ -2949,7 +2930,7 @@ FE_Nedelec::get_subface_interpolation_matrix( } } -#endif + // Since this is a vector valued element, // we cannot interpolate a scalar function. @@ -3811,7 +3792,7 @@ FE_Nedelec::interpolate (std::vector& local_dofs, // for the horizontal shape // functions. system_rhs = 0; - + for (unsigned int q_point = 0; q_point < n_face_points; ++q_point) { @@ -5447,6 +5428,9 @@ FE_Nedelec::memory_consumption () const } -template class FE_Nedelec; + +// explicit instantiations +#include "fe_nedelec.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_nedelec.inst.in b/deal.II/source/fe/fe_nedelec.inst.in new file mode 100644 index 0000000000..72a5fe8d82 --- /dev/null +++ b/deal.II/source/fe/fe_nedelec.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_Nedelec; + } + diff --git a/deal.II/source/fe/fe_nothing.cc b/deal.II/source/fe/fe_nothing.cc index 795e054e4c..debdf705d2 100644 --- a/deal.II/source/fe/fe_nothing.cc +++ b/deal.II/source/fe/fe_nothing.cc @@ -185,8 +185,8 @@ compare_for_face_domination (const FiniteElement & fe_other) const else return FiniteElementDomination::this_element_dominates; } - - + + template std::vector > FE_Nothing :: @@ -212,7 +212,7 @@ hp_line_dof_identities (const FiniteElement &/*fe_other*/) const return std::vector > (); } - + template std::vector > FE_Nothing :: @@ -225,7 +225,7 @@ hp_quad_dof_identities (const FiniteElement &/*fe_other*/) const return std::vector > (); } - + template bool FE_Nothing :: @@ -233,8 +233,8 @@ hp_constraints_are_implemented () const { return true; } - - + + template void FE_Nothing:: @@ -243,7 +243,7 @@ get_face_interpolation_matrix (const FiniteElement &/*source_fe*/, { // since this element has no face dofs, the // interpolation matrix is necessarily empty - + Assert (interpolation_matrix.m() == 0, ExcDimensionMismatch (interpolation_matrix.m(), 0)); @@ -271,9 +271,11 @@ get_subface_interpolation_matrix (const FiniteElement & /*source_fe*/, 0)); } - -template class FE_Nothing; + +// explicit instantiations +#include "fe_nothing.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_nothing.inst.in b/deal.II/source/fe/fe_nothing.inst.in new file mode 100644 index 0000000000..b2eb93ec76 --- /dev/null +++ b/deal.II/source/fe/fe_nothing.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_Nothing; + } + diff --git a/deal.II/source/fe/fe_poly.cc b/deal.II/source/fe/fe_poly.cc index d74db073ef..d4ae8dd54f 100644 --- a/deal.II/source/fe/fe_poly.cc +++ b/deal.II/source/fe/fe_poly.cc @@ -21,15 +21,6 @@ DEAL_II_NAMESPACE_OPEN -template class FE_Poly, deal_II_dimension>; -template class FE_Poly, deal_II_dimension>; -template class FE_Poly, deal_II_dimension>; +#include "fe_poly.inst" -#if deal_II_dimension != 3 - -template class FE_Poly, deal_II_dimension, deal_II_dimension+1>; -template class FE_Poly, deal_II_dimension, deal_II_dimension+1>; -//template class FE_Poly, deal_II_dimension, deal_II_dimension+1>; - -#endif DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_poly.inst.in b/deal.II/source/fe/fe_poly.inst.in new file mode 100644 index 0000000000..2a384acfa0 --- /dev/null +++ b/deal.II/source/fe/fe_poly.inst.in @@ -0,0 +1,29 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_Poly, deal_II_dimension>; + template class FE_Poly, deal_II_dimension>; + template class FE_Poly, deal_II_dimension>; + +#if deal_II_dimension != 3 + + template class FE_Poly, deal_II_dimension, deal_II_dimension+1>; + template class FE_Poly, deal_II_dimension, deal_II_dimension+1>; + //template class FE_Poly, deal_II_dimension, deal_II_dimension+1>; + +#endif + } + diff --git a/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/source/fe/fe_poly_tensor.cc index 7a313ba53d..0f58ad232d 100644 --- a/deal.II/source/fe/fe_poly_tensor.cc +++ b/deal.II/source/fe/fe_poly_tensor.cc @@ -35,7 +35,6 @@ namespace * algorithm, which determines the DoFs, which need this sign change for a * given cell. */ -#if deal_II_dimension == 1 void get_face_sign_change (const Triangulation<1>::cell_iterator &, const unsigned int , @@ -44,9 +43,9 @@ namespace // nothing to do in 1d std::fill (face_sign.begin (), face_sign.end (), 1.0); } -#endif -#if deal_II_dimension == 2 + + void get_face_sign_change (const Triangulation<2>::cell_iterator &cell, const unsigned int dofs_per_face, @@ -54,19 +53,19 @@ namespace { const unsigned int dim = 2; const unsigned int spacedim = 2; - + // Default is no sign // change. I.e. multiply by one. std::fill (face_sign.begin (), face_sign.end (), 1.0); - - for (unsigned int f = GeometryInfo::faces_per_cell / 2; + + for (unsigned int f = GeometryInfo::faces_per_cell / 2; f < GeometryInfo::faces_per_cell; ++f) { Triangulation::face_iterator face = cell->face (f); if (!face->at_boundary ()) { const unsigned int nn = cell->neighbor_face_no(f); - + if (nn < GeometryInfo::faces_per_cell / 2) for (unsigned int j = 0; j < dofs_per_face; ++j) { @@ -79,19 +78,17 @@ namespace } } } -#endif -#if deal_II_dimension == 3 + void get_face_sign_change (const Triangulation<3>::cell_iterator &/*cell*/, const unsigned int /*dofs_per_face*/, std::vector &face_sign) { std::fill (face_sign.begin (), face_sign.end (), 1.0); -//TODO: think about what it would take here +//TODO: think about what it would take here } -#endif } @@ -141,14 +138,14 @@ FE_PolyTensor::shape_value_component (const unsigned int i, { Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component < dim, ExcIndexRange (component, 0, dim)); - + if (cached_point != p || cached_values.size() == 0) { cached_point = p; cached_values.resize(poly_space.n()); poly_space.compute(p, cached_values, cached_grads, cached_grad_grads); } - + double s = 0; if (inverse_node_matrix.n_cols() == 0) return cached_values[i][component]; @@ -180,21 +177,21 @@ FE_PolyTensor::shape_grad_component (const unsigned int i, { Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component < dim, ExcIndexRange (component, 0, dim)); - + if (cached_point != p || cached_grads.size() == 0) { cached_point = p; cached_grads.resize(poly_space.n()); poly_space.compute(p, cached_values, cached_grads, cached_grad_grads); } - + Tensor<1,dim> s; if (inverse_node_matrix.n_cols() == 0) return cached_grads[i][component]; else for (unsigned int j=0;j::shape_grad_grad_component (const unsigned int { Assert (idofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell)); Assert (component < dim, ExcIndexRange (component, 0, dim)); - + if (cached_point != p || cached_grad_grads.size() == 0) { cached_point = p; cached_grad_grads.resize(poly_space.n()); poly_space.compute(p, cached_values, cached_grads, cached_grad_grads); } - + Tensor<2,dim> s; if (inverse_node_matrix.n_cols() == 0) return cached_grad_grads[i][component]; else for (unsigned int j=0;j::get_data ( const UpdateFlags flags(data->update_flags); const unsigned int n_q_points = quadrature.size(); - + // some scratch arrays std::vector > values(0); std::vector > grads(0); @@ -295,7 +292,7 @@ FE_PolyTensor::get_data ( // that if (flags & update_hessians) { -// grad_grads.resize (this->dofs_per_cell); +// grad_grads.resize (this->dofs_per_cell); data->initialize_2nd (this, mapping, quadrature); } @@ -310,7 +307,7 @@ FE_PolyTensor::get_data ( { poly_space.compute(quadrature.point(k), values, grads, grad_grads); - + if (flags & update_values) { if (inverse_node_matrix.n_cols() == 0) @@ -325,7 +322,7 @@ FE_PolyTensor::get_data ( data->shape_values[i][k] = add_values; } } - + if (flags & update_gradients) { if (inverse_node_matrix.n_cols() == 0) @@ -340,7 +337,7 @@ FE_PolyTensor::get_data ( data->shape_grads[i][k] = add_grads; } } - + } return data; } @@ -373,7 +370,7 @@ FE_PolyTensor::fill_fe_values ( const unsigned int n_q_points = quadrature.size(); const UpdateFlags flags(fe_data.current_update_flags()); - + Assert(!(flags & update_values) || fe_data.shape_values.size() == this->dofs_per_cell, ExcDimensionMismatch(fe_data.shape_values.size(), this->dofs_per_cell)); Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_q_points, @@ -394,11 +391,11 @@ FE_PolyTensor::fill_fe_values ( if (mapping_type == mapping_piola) if (cell_similarity == CellSimilarity::translation) cell_similarity = CellSimilarity::none; - + for (unsigned int i=0; idofs_per_cell; ++i) { const unsigned int first = data.shape_function_to_row_table[i]; - + if (flags & update_values && cell_similarity != CellSimilarity::translation) switch (mapping_type) { @@ -409,7 +406,7 @@ FE_PolyTensor::fill_fe_values ( data.shape_values(first+d,k) = fe_data.shape_values[i][k][d]; break; } - + case mapping_covariant: case mapping_contravariant: { @@ -434,24 +431,24 @@ FE_PolyTensor::fill_fe_values ( = sign_change[i] * shape_values[k][d]; break; } - + case mapping_nedelec: { std::vector > shape_values (n_q_points); mapping.transform (fe_data.shape_values[i], shape_values, mapping_data, mapping_covariant); - + for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) data.shape_values(first+d,k) = shape_values[k][d]; - + break; } default: Assert(false, ExcNotImplemented()); } - + if (flags & update_gradients && cell_similarity != CellSimilarity::translation) { std::vector > shape_grads1 (n_q_points); @@ -472,36 +469,36 @@ FE_PolyTensor::fill_fe_values ( { mapping.transform(fe_data.shape_grads[i], shape_grads1, mapping_data, mapping_covariant_gradient); - + for (unsigned int k=0; k::fill_fe_values ( // (J^-1) from the right, // we have to trick a // little in between - + // do first transformation mapping.transform (fe_data.shape_grads[i], shape_grads1, mapping_data, mapping_covariant); - // transpose matrix + // transpose matrix for (unsigned int k = 0; k < n_q_points; ++k) shape_grads2[k] = transpose (shape_grads1[k]); - // do second transformation + // do second transformation mapping.transform (shape_grads2, shape_grads1, mapping_data, mapping_covariant); - // transpose back + // transpose back for (unsigned int k = 0; k < n_q_points; ++k) shape_grads2[k] = transpose (shape_grads1[k]); - + for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_gradients[first + d][k] = shape_grads2[k][d]; + data.shape_gradients[first + d][k] = shape_grads2[k][d]; // then copy over to target: break; } - + default: Assert(false, ExcNotImplemented()); } } } - + if (flags & update_hessians && cell_similarity != CellSimilarity::translation) this->compute_2nd (mapping, cell, typename QProjector::DataSetDescriptor().cell(), @@ -579,7 +576,7 @@ FE_PolyTensor::fill_fe_face_values ( // offset determines which data set // to take (all data sets for all // faces are stored contiguously) - + const typename QProjector::DataSetDescriptor dsd; const typename QProjector::DataSetDescriptor offset = dsd.face (face, @@ -587,7 +584,7 @@ FE_PolyTensor::fill_fe_face_values ( cell->face_flip(face), cell->face_rotation(face), n_q_points); - + const UpdateFlags flags(fe_data.update_once | fe_data.update_each); //TODO: Size assertions @@ -599,11 +596,11 @@ FE_PolyTensor::fill_fe_face_values ( // on the neighborhood between two faces. std::vector sign_change (this->dofs_per_cell, 1.0); get_face_sign_change (cell, this->dofs_per_face, sign_change); - + for (unsigned int i=0; idofs_per_cell; ++i) { const unsigned int first = data.shape_function_to_row_table[i]; - + if (flags & update_values) { switch (mapping_type) @@ -615,7 +612,7 @@ FE_PolyTensor::fill_fe_face_values ( data.shape_values(first+d,k) = fe_data.shape_values[i][k+offset][d]; break; } - + case mapping_covariant: case mapping_contravariant: { @@ -624,11 +621,11 @@ FE_PolyTensor::fill_fe_face_values ( std::vector > shape_values (n_q_points); mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), shape_values, mapping_data, mapping_type); - + for (unsigned int k=0; k::fill_fe_face_values ( = sign_change[i] * shape_values[k][d]; break; } - + case mapping_nedelec: { std::vector > shape_values (n_q_points); mapping.transform (make_slice (fe_data.shape_values[i], offset, n_q_points), shape_values, mapping_data, mapping_covariant); - + for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) data.shape_values(first+d,k) = shape_values[k][d]; - + break; } - + default: Assert(false, ExcNotImplemented()); - } + } } - + if (flags & update_gradients) { std::vector > shape_grads1 (n_q_points); @@ -676,12 +673,12 @@ FE_PolyTensor::fill_fe_face_values ( data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - + case mapping_covariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant_gradient); - + for (unsigned int k=0; k::fill_fe_face_values ( for (unsigned int k=0; k::fill_fe_face_values ( // (J^-1) from the right, // we have to trick a // little in between - + // do first transformation mapping.transform (make_slice (fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant); - // transpose matrix + // transpose matrix for (unsigned int k = 0; k < n_q_points; ++k) shape_grads2[k] = transpose (shape_grads1[k]); - // do second transformation + // do second transformation mapping.transform (shape_grads2, shape_grads1, mapping_data, mapping_covariant); - // transpose back + // transpose back for (unsigned int k = 0; k < n_q_points; ++k) shape_grads2[k] = transpose (shape_grads1[k]); - + for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_gradients[first + d][k] = shape_grads2[k][d]; + data.shape_gradients[first + d][k] = shape_grads2[k][d]; // then copy over to target: break; } - + default: Assert(false, ExcNotImplemented()); } } } - + if (flags & update_hessians) this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data); } @@ -810,11 +807,11 @@ FE_PolyTensor::fill_fe_subface_values ( // on the neighborhood between two faces. std::vector sign_change (this->dofs_per_cell, 1.0); get_face_sign_change (cell, this->dofs_per_face, sign_change); - + for (unsigned int i=0; idofs_per_cell; ++i) { const unsigned int first = data.shape_function_to_row_table[i]; - + if (flags & update_values) { switch (mapping_type) @@ -826,7 +823,7 @@ FE_PolyTensor::fill_fe_subface_values ( data.shape_values(first+d,k) = fe_data.shape_values[i][k+offset][d]; break; } - + case mapping_covariant: case mapping_contravariant: { @@ -835,11 +832,11 @@ FE_PolyTensor::fill_fe_subface_values ( std::vector > shape_values (n_q_points); mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points), shape_values, mapping_data, mapping_type); - + for (unsigned int k=0; k::fill_fe_subface_values ( std::vector > shape_values (n_q_points); mapping.transform (make_slice (fe_data.shape_values[i], offset, n_q_points), shape_values, mapping_data, mapping_covariant); - + for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) data.shape_values(first+d,k) = shape_values[k][d]; - + break; } default: Assert(false, ExcNotImplemented()); } } - + if (flags & update_gradients) { std::vector > shape_grads1 (n_q_points); @@ -885,7 +882,7 @@ FE_PolyTensor::fill_fe_subface_values ( data.shape_gradients[first+d][k] = shape_grads1[k][d]; break; } - + case mapping_covariant: { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), @@ -901,11 +898,11 @@ FE_PolyTensor::fill_fe_subface_values ( { mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_contravariant_gradient); - + for (unsigned int k=0; k::fill_fe_subface_values ( mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_piola); - + for (unsigned int k=0; k::fill_fe_subface_values ( // (J^-1) from the right, // we have to trick a // little in between - + // do first transformation mapping.transform (make_slice (fe_data.shape_grads[i], offset, n_q_points), shape_grads1, mapping_data, mapping_covariant); - // transpose matrix + // transpose matrix for (unsigned int k = 0; k < n_q_points; ++k) shape_grads2[k] = transpose (shape_grads1[k]); - // do second transformation + // do second transformation mapping.transform (shape_grads2, shape_grads1, mapping_data, mapping_covariant); - // transpose back + // transpose back for (unsigned int k = 0; k < n_q_points; ++k) shape_grads2[k] = transpose (shape_grads1[k]); - + for (unsigned int k = 0; k < n_q_points; ++k) for (unsigned int d = 0; d < dim; ++d) - data.shape_gradients[first + d][k] = shape_grads2[k][d]; + data.shape_gradients[first + d][k] = shape_grads2[k][d]; // then copy over to target: break; } - + default: Assert(false, ExcNotImplemented()); } } } - + if (flags & update_hessians) this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data); } @@ -1007,7 +1004,7 @@ UpdateFlags FE_PolyTensor::update_once (const UpdateFlags flags) const { const bool values_once = (mapping_type == mapping_none); - + UpdateFlags out = update_default; if (values_once && (flags & update_values)) @@ -1029,26 +1026,26 @@ FE_PolyTensor::update_each (const UpdateFlags flags) const { if (flags & update_values) out |= update_values | update_piola; - + if (flags & update_gradients) out |= update_gradients | update_piola | update_covariant_transformation; - + if (flags & update_hessians) out |= update_hessians | update_piola | update_covariant_transformation; - + break; } - + case mapping_nedelec: { if (flags & update_values) out |= update_values | update_covariant_transformation; - + if (flags & update_gradients) out |= update_gradients | update_covariant_transformation; - + if (flags & update_hessians) out |= update_hessians | update_covariant_transformation; - + break; } @@ -1057,16 +1054,14 @@ FE_PolyTensor::update_each (const UpdateFlags flags) const Assert (false, ExcNotImplemented()); } } - + return out; } -template class FE_PolyTensor,deal_II_dimension>; -template class FE_PolyTensor,deal_II_dimension>; -template class FE_PolyTensor,deal_II_dimension>; -template class FE_PolyTensor, deal_II_dimension>; +// explicit instantiations +#include "fe_poly_tensor.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_poly_tensor.inst.in b/deal.II/source/fe/fe_poly_tensor.inst.in new file mode 100644 index 0000000000..945f4670d7 --- /dev/null +++ b/deal.II/source/fe/fe_poly_tensor.inst.in @@ -0,0 +1,22 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_PolyTensor,deal_II_dimension>; + template class FE_PolyTensor,deal_II_dimension>; + template class FE_PolyTensor,deal_II_dimension>; + template class FE_PolyTensor, deal_II_dimension>; + } + diff --git a/deal.II/source/fe/fe_q.cc b/deal.II/source/fe/fe_q.cc index e924afe60c..c1d4f0ddfa 100644 --- a/deal.II/source/fe/fe_q.cc +++ b/deal.II/source/fe/fe_q.cc @@ -627,7 +627,7 @@ FE_Q::FE_Q (const Quadrature<1> &points) // embedding operators FETools::compute_embedding_matrices (*this, this->prolongation); - // Fill restriction matrices + // Fill restriction matrices initialize_restriction(); initialize_quad_dof_index_permutation(); @@ -802,7 +802,6 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, } -#if deal_II_dimension == 1 template <> void @@ -856,10 +855,7 @@ get_subface_interpolation_matrix (const FiniteElement<1,2> &/*x_source_fe*/, ExcInterpolationNotImplemented ()); } -#endif - -#if deal_II_dimension > 1 template void @@ -1076,7 +1072,6 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe } } -#endif template @@ -1329,7 +1324,6 @@ void FE_Q::initialize_unit_support_points (const Quadrature<1> &po } -#if deal_II_dimension == 1 template <> void FE_Q<1>::initialize_unit_face_support_points () @@ -1355,7 +1349,6 @@ void FE_Q<1,2>::initialize_unit_face_support_points (const Quadrature<1> &/*poin // no faces in 1d, so nothing to do } -#endif template @@ -1437,7 +1430,6 @@ FE_Q::initialize_quad_dof_index_permutation () } -#if deal_II_dimension == 3 template <> void @@ -1505,7 +1497,6 @@ FE_Q<3>::initialize_quad_dof_index_permutation () this->adjust_line_dof_index_for_line_orientation_table[i]=this->dofs_per_line-1-i - i; } -#endif @@ -1531,7 +1522,6 @@ FE_Q::face_lexicographic_to_hierarchic_numbering (const unsigned i } -#if deal_II_dimension == 1 template <> std::vector @@ -1547,7 +1537,7 @@ FE_Q<1,2>::face_lexicographic_to_hierarchic_numbering (const unsigned int) return std::vector(); } -#endif + template void @@ -1811,7 +1801,7 @@ FE_Q::initialize_restriction () // value per row. Still, the // values should sum up to 1 double sum_check = 0; - for (unsigned int child_dof = 0; child_dofdofs_per_cell; + for (unsigned int child_dof = 0; child_dofdofs_per_cell; ++child_dof) { const double val @@ -1823,7 +1813,7 @@ FE_Q::initialize_restriction () sum_check += val; } - Assert (std::fabs(sum_check-1) < this->degree*eps, + Assert (std::fabs(sum_check-1) < this->degree*eps, ExcInternalError()); } } @@ -1831,11 +1821,11 @@ FE_Q::initialize_restriction () } + //--------------------------------------------------------------------------- // Data field initialization //--------------------------------------------------------------------------- -#if deal_II_dimension == 1 template <> bool @@ -1882,10 +1872,6 @@ FE_Q<1,2>::has_support_on_face (const unsigned int shape_index, - -#else - - template bool FE_Q::has_support_on_face (const unsigned int shape_index, @@ -1992,7 +1978,6 @@ FE_Q::has_support_on_face (const unsigned int shape_index, return false; } -#endif template @@ -2005,12 +1990,7 @@ FE_Q::memory_consumption () const -template class FE_Q; - -#if deal_II_dimension !=3 -template class FE_Q; -#endif - - +// explicit instantiations +#include "fe_q.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_q.inst.in b/deal.II/source/fe/fe_q.inst.in new file mode 100644 index 0000000000..d72cebf3a4 --- /dev/null +++ b/deal.II/source/fe/fe_q.inst.in @@ -0,0 +1,23 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_Q; + +#if deal_II_dimension !=3 + template class FE_Q; +#endif + } + diff --git a/deal.II/source/fe/fe_q_hierarchical.cc b/deal.II/source/fe/fe_q_hierarchical.cc index f3b21f9134..64cebd601b 100644 --- a/deal.II/source/fe/fe_q_hierarchical.cc +++ b/deal.II/source/fe/fe_q_hierarchical.cc @@ -24,7 +24,7 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace { inline std::vector @@ -53,20 +53,20 @@ FE_Q_Hierarchical::FE_Q_Hierarchical (const unsigned int degree) { this->poly_space.set_numbering( hierarchic_to_fe_q_hierarchical_numbering(*this)); - - // The matrix @p{dofs_cell} contains the - // values of the linear functionals of - // the master 1d cell applied to the + + // The matrix @p{dofs_cell} contains the + // values of the linear functionals of + // the master 1d cell applied to the // shape functions of the two 1d subcells. // The matrix @p{dofs_subcell} constains - // the values of the linear functionals - // on each 1d subcell applied to the - // shape functions on the master 1d - // subcell. - // We use @p{dofs_cell} and - // @p{dofs_subcell} to compute the - // @p{prolongation}, @p{restriction} and - // @p{interface_constraints} matrices + // the values of the linear functionals + // on each 1d subcell applied to the + // shape functions on the master 1d + // subcell. + // We use @p{dofs_cell} and + // @p{dofs_subcell} to compute the + // @p{prolongation}, @p{restriction} and + // @p{interface_constraints} matrices // for all dimensions. std::vector > dofs_cell (GeometryInfo<1>::max_children_per_cell, @@ -105,7 +105,7 @@ FE_Q_Hierarchical::get_name () const // this function returns, so they // have to be kept in synch - std::ostringstream namebuf; + std::ostringstream namebuf; namebuf << "FE_Q_Hierarchical<" << dim << ">(" << this->degree << ")"; return namebuf.str(); @@ -183,12 +183,12 @@ FE_Q_Hierarchical::build_dofs_cell (std::vector > &dofs_ // upper diagonal block if ((j<=1) && (k<=1)) { - if (((c==0) && (j==0) && (k==0)) || + if (((c==0) && (j==0) && (k==0)) || ((c==1) && (j==1) && (k==1))) dofs_cell[c](j,k) = 1.; - else + else dofs_cell[c](j,k) = 0.; - + if (((c==0) && (j==1)) || ((c==1) && (j==0))) dofs_subcell[c](j,k) = .5; else if (((c==0) && (k==0)) || ((c==1) && (k==1))) @@ -211,7 +211,7 @@ FE_Q_Hierarchical::build_dofs_cell (std::vector > &dofs_ factor *= ((double) (k-i+1))/((double) i); if (c==0) { - dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? + dofs_subcell[c](j,k) = ((k+j) % 2 == 0) ? std::pow(.5,static_cast(k))*factor : -std::pow(.5,static_cast(k))*factor; dofs_cell[c](j,k) = std::pow(2.,static_cast(j))*factor; @@ -219,7 +219,7 @@ FE_Q_Hierarchical::build_dofs_cell (std::vector > &dofs_ else { dofs_subcell[c](j,k) = std::pow(.5,static_cast(k))*factor; - dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? + dofs_cell[c](j,k) = ((k+j) % 2 == 0) ? std::pow(2.,static_cast(j))*factor : -std::pow(2.,static_cast(j))*factor; } @@ -247,58 +247,58 @@ initialize_constraints (const std::vector > &dofs_subcell) // no constraints in 1d break; } - + case 2: { // vertex node for (unsigned int i=0; iinterface_constraints(0,i) = dofs_subcell[0](1,i); + this->interface_constraints(0,i) = dofs_subcell[0](1,i); // edge nodes for (unsigned int c=0; c::max_children_per_cell; ++c) for (unsigned int i=0; iinterface_constraints(1 + c*(degree-1) + j - 2,i) = + this->interface_constraints(1 + c*(degree-1) + j - 2,i) = dofs_subcell[c](j,i); break; } - + case 3: { for (unsigned int i=0; iinterface_constraints(0,face_renumber[i]) = - dofs_subcell[0](1,i % dofs_1d) * + // center vertex node + this->interface_constraints(0,face_renumber[i]) = + dofs_subcell[0](1,i % dofs_1d) * dofs_subcell[0](1,(i - (i % dofs_1d)) / dofs_1d); // boundary vertex nodes - this->interface_constraints(1,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * + this->interface_constraints(1,face_renumber[i]) = + dofs_subcell[0](0, i % dofs_1d) * dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(2,face_renumber[i]) = - dofs_subcell[1](1, i % dofs_1d) * + this->interface_constraints(2,face_renumber[i]) = + dofs_subcell[1](1, i % dofs_1d) * dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(3,face_renumber[i]) = - dofs_subcell[0](1, i % dofs_1d) * + this->interface_constraints(3,face_renumber[i]) = + dofs_subcell[0](1, i % dofs_1d) * dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(4,face_renumber[i]) = - dofs_subcell[1](0, i % dofs_1d) * + this->interface_constraints(4,face_renumber[i]) = + dofs_subcell[1](0, i % dofs_1d) * dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); - + // interior edges for (unsigned int j=0; j<(degree-1); j++) { - this->interface_constraints(5 + j,face_renumber[i]) = - dofs_subcell[0](1, i % dofs_1d) * + this->interface_constraints(5 + j,face_renumber[i]) = + dofs_subcell[0](1, i % dofs_1d) * dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = - dofs_subcell[0](1,i % dofs_1d) * + this->interface_constraints(5 + (degree-1) + j,face_renumber[i]) = + dofs_subcell[0](1,i % dofs_1d) * dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j,i % dofs_1d) * + this->interface_constraints(5 + 2*(degree-1) + j,face_renumber[i]) = + dofs_subcell[0](2 + j,i % dofs_1d) * dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d); - this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * + this->interface_constraints(5 + 3*(degree-1) + j,face_renumber[i]) = + dofs_subcell[1](2 + j, i % dofs_1d) * dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d); } @@ -307,31 +307,31 @@ initialize_constraints (const std::vector > &dofs_subcell) { // left edge this->interface_constraints(5 + 4*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[0](0, i % dofs_1d) * dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + (degree-1) + j,face_renumber[i]) = - dofs_subcell[0](0, i % dofs_1d) * + dofs_subcell[0](0, i % dofs_1d) * dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); // right edge this->interface_constraints(5 + 4*(degree-1) + 2*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](1, i % dofs_1d) * + dofs_subcell[1](1, i % dofs_1d) * dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + 3*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](1, i % dofs_1d) * + dofs_subcell[1](1, i % dofs_1d) * dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d); - // bottom edge + // bottom edge this->interface_constraints(5 + 4*(degree-1) + 4*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + j, i % dofs_1d) * dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + 5*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + j, i % dofs_1d) * dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d); // top edge this->interface_constraints(5 + 4*(degree-1) + 6*(degree-1) + j,face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + j, i % dofs_1d) * dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); this->interface_constraints(5 + 4*(degree-1) + 7*(degree-1) + j,face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + j, i % dofs_1d) * dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d); } @@ -341,25 +341,25 @@ initialize_constraints (const std::vector > &dofs_subcell) { // subcell 0 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1),face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + j, i % dofs_1d) * dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); // subcell 1 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + (degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + j, i % dofs_1d) * dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d); // subcell 2 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 2*(degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[0](2 + j, i % dofs_1d) * + dofs_subcell[0](2 + j, i % dofs_1d) * dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); // subcell 3 this->interface_constraints(5 + 12*(degree-1) + j + k*(degree-1) + 3*(degree-1)*(degree-1),face_renumber[i]) = - dofs_subcell[1](2 + j, i % dofs_1d) * + dofs_subcell[1](2 + j, i % dofs_1d) * dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d); } } break; } - + default: Assert (false, ExcNotImplemented()); } @@ -378,7 +378,7 @@ initialize_embedding_and_restriction (const std::vector > &do const unsigned int dofs_1d = 2*this->dofs_per_vertex + this->dofs_per_line; const std::vector &renumber= this->poly_space.get_numbering(); - + for (unsigned int c=0; c::max_children_per_cell; ++c) { this->prolongation[iso][c].reinit (this->dofs_per_cell, this->dofs_per_cell); @@ -399,7 +399,7 @@ initialize_embedding_and_restriction (const std::vector > &do // for higher dimensions, things // are a little more tricky: - + // j loops over dofs in the // subcell. These are the rows in // the embedding matrix. @@ -420,13 +420,13 @@ initialize_embedding_and_restriction (const std::vector > &do // bottom/top line: 0/1 const unsigned int c1 = c/2; - this->prolongation[iso][c](j,i) = + this->prolongation[iso][c](j,i) = dofs_subcell[c0](renumber[j] % dofs_1d, renumber[i] % dofs_1d) * dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d, (renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d); - this->restriction[iso][c](j,i) = + this->restriction[iso][c](j,i) = dofs_cell[c0](renumber[j] % dofs_1d, renumber[i] % dofs_1d) * dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d, @@ -434,7 +434,7 @@ initialize_embedding_and_restriction (const std::vector > &do } break; } - + case 3: { for (unsigned int c=0; c::max_children_per_cell; ++c) @@ -446,7 +446,7 @@ initialize_embedding_and_restriction (const std::vector > &do // bottom/top face: 0/1 const unsigned int c2 = c/4; - this->prolongation[iso][c](j,i) = + this->prolongation[iso][c](j,i) = dofs_subcell[c0](renumber[j] % dofs_1d, renumber[i] % dofs_1d) * dofs_subcell[c1](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) % dofs_1d, @@ -454,7 +454,7 @@ initialize_embedding_and_restriction (const std::vector > &do dofs_subcell[c2](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d - (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d, ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d - (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d ) % dofs_1d)) / dofs_1d); - this->restriction[iso][c](j,i) = + this->restriction[iso][c](j,i) = dofs_cell[c0](renumber[j] % dofs_1d, renumber[i] % dofs_1d) * dofs_cell[c1](((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) % dofs_1d, @@ -479,12 +479,12 @@ void FE_Q_Hierarchical::initialize_unit_support_points () unsigned int n = this->degree+1; for (unsigned int i=1; idegree+1; - + this->unit_support_points.resize(n); const std::vector &index_map_inverse= this->poly_space.get_numbering_inverse(); - + Point p; // the method of numbering allows // each dof to be associated with a @@ -544,7 +544,6 @@ void FE_Q_Hierarchical::initialize_unit_support_points () } -#if deal_II_dimension == 1 template <> void FE_Q_Hierarchical<1>::initialize_unit_face_support_points () @@ -576,11 +575,8 @@ get_subface_interpolation_matrix (const FiniteElement<1> &/*x_source_fe*/, ExcInterpolationNotImplemented ()); } -#endif -#if deal_II_dimension > 1 - template void FE_Q_Hierarchical:: @@ -625,25 +621,25 @@ get_face_interpolation_matrix (const FiniteElement &x_source_fe, typename FEL:: ExcInterpolationNotImplemented ()); interpolation_matrix = 0; - + switch (dim) { case 2: { for (unsigned int i = 0; i < this->dofs_per_face; ++i) interpolation_matrix (i, i) = 1; - + break; } - + case 3: { for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i) interpolation_matrix (i, i) = 1; - + for (unsigned int i = 0; i < this->degree - 1; ++i) { for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j) interpolation_matrix ( i + j * (x_source_fe.degree - 1) + GeometryInfo<3>::vertices_per_face, i + j * (this->degree - 1) + GeometryInfo<3>::vertices_per_face) = 1; - + for (unsigned int j = 0; j < this->degree - 1; ++j) interpolation_matrix ( (i + GeometryInfo<3>::lines_per_face) * (x_source_fe.degree - 1) + j @@ -701,7 +697,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, Assert (this->dofs_per_face <= source_fe.dofs_per_face, typename FEL:: ExcInterpolationNotImplemented ()); - + switch (dim) { case 2: @@ -713,68 +709,68 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (0, 0) = 1.0; interpolation_matrix (1, 0) = 0.5; interpolation_matrix (1, 1) = 0.5; - + for (unsigned int dof = 2; dof < this->dofs_per_face;) { interpolation_matrix (1, dof) = -1.0; dof = dof + 2; } - + int factorial_i = 1; int factorial_ij; int factorial_j; - + for (int i = 2; i < (int) this->dofs_per_face; ++i) { interpolation_matrix (i, i) = std::pow (0.5, i); factorial_i *= i; factorial_j = factorial_i; factorial_ij = 1; - + for (int j = i + 1; j < (int) this->dofs_per_face; ++j) { factorial_ij *= j - i; factorial_j *= j; - + if ((i + j) & 1) interpolation_matrix (i, j) = -1.0 * std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij); - + else interpolation_matrix (i, j) = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij); } } - + break; } - + case 1: { interpolation_matrix (0, 0) = 0.5; interpolation_matrix (0, 1) = 0.5; - + for (unsigned int dof = 2; dof < this->dofs_per_face;) { interpolation_matrix (0, dof) = -1.0; dof = dof + 2; } - + interpolation_matrix (1, 1) = 1.0; - + int factorial_i = 1; int factorial_ij; int factorial_j; - + for (int i = 2; i < (int) this->dofs_per_face; ++i) { interpolation_matrix (i, i) = std::pow (0.5, i); factorial_i *= i; factorial_j = factorial_i; factorial_ij = 1; - + for (int j = i + 1; j < (int) this->dofs_per_face; ++j) { factorial_ij *= j - i; @@ -785,10 +781,10 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, } } } - + break; } - + case 3: { switch (subface) @@ -800,33 +796,33 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (1, 1) = 0.5; interpolation_matrix (2, 0) = 0.5; interpolation_matrix (2, 2) = 0.5; - + for (unsigned int i = 0; i < this->degree - 1;) { for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) interpolation_matrix (3, i + line * (this->degree - 1) + 4) = -0.5; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (3, i + (j + 4) * this->degree - j) = 1.0; j = j + 2; } - + interpolation_matrix (1, i + 2 * (this->degree + 1)) = -1.0; interpolation_matrix (2, i + 4) = -1.0; i = i + 2; } - + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) interpolation_matrix (3, vertex) = 0.25; - + int factorial_i = 1; int factorial_ij; int factorial_j; int factorial_k; int factorial_kl; int factorial_l; - + for (int i = 2; i <= (int) this->degree; ++i) { double tmp = std::pow (0.5, i); @@ -838,78 +834,78 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; tmp *= -2.0; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (i + source_fe.degree + 1, (i + 2) * this->degree + j + 2 - i) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, i + (j + 4) * this->degree - j - 2) = tmp; j = j + 2; } - + factorial_k = 1; - + for (int j = 2; j <= (int) this->degree; ++j) { interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); factorial_k *= j; factorial_kl = 1; factorial_l = factorial_k; - + for (int k = j + 1; k < (int) this->degree; ++k) { factorial_kl *= k - j; factorial_l *= k; - + if ((j + k) & 1) interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = -1.0 * std::pow (0.5, i + k) * factorial_l / (factorial_k * factorial_kl); - + else interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } - + factorial_i *= i; factorial_j = factorial_i; factorial_ij = 1; - + for (int j = i + 1; j <= (int) this->degree; ++j) { factorial_ij *= j - i; factorial_j *= j; - + if ((i + j) & 1) { tmp = -1.0 * std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix (i + 2, j + 2) = tmp; interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; factorial_k = 1; - + for (int k = 2; k <= (int) this->degree; ++k) { interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); factorial_k *= k; factorial_l = factorial_k; factorial_kl = 1; - + for (int l = k + 1; l <= (int) this->degree; ++l) { factorial_kl *= l - k; factorial_l *= l; - + if ((k + l) & 1) interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = -1.0 * tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); - + else interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); } } - + tmp *= 0.5; interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp; interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; tmp *= -2.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + source_fe.degree + 1, (j + 2) * this->degree + k + 2 - j) = tmp; @@ -923,34 +919,34 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + 2, j + 2) = tmp; interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; factorial_k = 1; - + for (int k = 2; k <= (int) this->degree; ++k) { interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); factorial_k *= k; factorial_l = factorial_k; factorial_kl = 1; - + for (int l = k + 1; l <= (int) this->degree; ++l) { factorial_kl *= l - k; factorial_l *= l; - + if ((k + l) & 1) interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = -1.0 * tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); - + else interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); } } - + tmp *= 0.5; interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp; interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; tmp *= -2.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + source_fe.degree + 1, (j + 2) * this->degree + k + 2 - j) = tmp; @@ -960,10 +956,10 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, } } } - + break; } - + case 1: { interpolation_matrix (0, 0) = 0.5; @@ -971,33 +967,33 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (1, 1) = 1.0; interpolation_matrix (3, 1) = 0.5; interpolation_matrix (3, 3) = 0.5; - + for (unsigned int i = 0; i < this->degree - 1;) { for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) interpolation_matrix (2, i + line * (this->degree - 1) + 4) = -0.5; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (2, i + (j + 4) * this->degree - j) = 1.0; j = j + 2; } - + interpolation_matrix (0, i + 2 * (this->degree + 1)) = -1.0; interpolation_matrix (3, i + this->degree + 3) = -1.0; i = i + 2; } - + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) interpolation_matrix (2, vertex) = 0.25; - + int factorial_i = 1; int factorial_ij; int factorial_j; int factorial_k; int factorial_kl; int factorial_l; - + for (int i = 2; i <= (int) this->degree; ++i) { double tmp = std::pow (0.5, i + 1); @@ -1006,21 +1002,21 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + 3 * source_fe.degree - 1, i + 2 * this->degree) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; tmp *= -2.0; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (i + 2, j + (i + 2) * this->degree + 2 - i) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, i + (j + 4) * this->degree - j - 2) = tmp; j = j + 2; } - + tmp *= - 1.0; interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp; interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; factorial_i *= i; factorial_j = factorial_i; factorial_ij = 1; - + for (int j = i + 1; j <= (int) this->degree; ++j) { factorial_ij *= j - i; @@ -1028,80 +1024,80 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, tmp = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; factorial_k = 1; - + for (int k = 2; k <= (int) this->degree; ++k) { interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); factorial_k *= k; factorial_l = factorial_k; factorial_kl = 1; - + for (int l = k + 1; l <= (int) this->degree; ++l) { factorial_kl *= l - k; factorial_l *= l; - + if ((k + l) & 1) interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = -1.0 * tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); - + else interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); } } - + tmp *= -1.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + 3 * source_fe.degree - 1, j + (k + 4) * this->degree - k - 2) = tmp; k = k + 2; } - + tmp *= -0.5; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 2 * this->degree) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; - + if ((i + j) & 1) tmp *= -1.0; - + interpolation_matrix (i + 2, j + 2) = tmp; interpolation_matrix (i + 2, j + this->degree + 1) = tmp; interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = 2.0 * tmp; tmp *= -2.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + 2, k + (j + 2) * this->degree + 2 - j) = tmp; k = k + 2; } } - + factorial_k = 1; - + for (int j = 2; j <= (int) this->degree; ++j) { interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); factorial_k *= j; factorial_l = factorial_k; factorial_kl = 1; - + for (int k = j + 1; k <= (int) this->degree; ++k) { factorial_kl *= k - j; factorial_l *= k; - + if ((j + k) & 1) interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = -1.0 * std::pow (0.5, i + k) * factorial_l / (factorial_k * factorial_kl); - + else interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } } - + break; } - + case 2: { interpolation_matrix (0, 0) = 0.5; @@ -1109,33 +1105,33 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (2, 2) = 1.0; interpolation_matrix (3, 2) = 0.5; interpolation_matrix (3, 3) = 0.5; - + for (unsigned int i = 0; i < this->degree - 1;) { for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) interpolation_matrix (1, i + line * (this->degree - 1) + 4) = -0.5; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (1, i + (j + 4) * this->degree - j) = 1.0; j = j + 2; } - + interpolation_matrix (0, i + 4) = -1.0; interpolation_matrix (3, i + 3 * this->degree + 1) = -1.0; i = i + 2; } - + for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) interpolation_matrix (1, vertex) = 0.25; - + int factorial_i = 1; int factorial_ij; int factorial_j; int factorial_k; int factorial_kl; int factorial_l; - + for (int i = 2; i <= (int) this->degree; ++i) { double tmp = std::pow (0.5, i); @@ -1147,23 +1143,23 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; interpolation_matrix (i + 2 * source_fe.degree, i + 3 * this->degree - 1) = tmp; tmp *= -2.0; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (i + source_fe.degree + 1, j + (i + 2) * this->degree + 2 - i) = tmp; interpolation_matrix (i + 2 * source_fe.degree, i + (j + 4) * this->degree - j - 2) = tmp; j = j + 2; } - + factorial_k = 1; - + for (int j = 2; j <= (int) this->degree; ++j) { interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); factorial_k *= j; factorial_l = factorial_k; factorial_kl = 1; - + for (int k = j + 1; k <= (int) this->degree; ++k) { factorial_kl *= k - j; @@ -1171,11 +1167,11 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } - + factorial_i *= i; factorial_j = factorial_i; factorial_ij = 1; - + for (int j = i + 1; j <= (int) this->degree; ++j) { factorial_ij *= j - i; @@ -1183,33 +1179,33 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, tmp = std::pow (0.5, j) * factorial_j / (factorial_i * factorial_ij); interpolation_matrix (i + 2, j + 2) = tmp; tmp *= -1.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + source_fe.degree + 1, k + (j + 2) * this->degree + 2 - j) = tmp; k = k + 2; } - + tmp *= -0.5; interpolation_matrix (i + source_fe.degree + 1, j + 2) = tmp; interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; - + if ((i + j) & 1) tmp *= -1.0; - + interpolation_matrix (i + 2 * source_fe.degree, j + 2 * this->degree) = tmp; interpolation_matrix (i + 2 * source_fe.degree, j + 3 * this->degree - 1) = tmp; tmp *= 2.0; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; factorial_k = 1; - + for (int k = 2; k <= (int) this->degree; ++k) { interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); factorial_k *= k; factorial_l = factorial_k; factorial_kl = 1; - + for (int l = k + 1; l <= (int) this->degree; ++l) { factorial_kl *= l - k; @@ -1217,9 +1213,9 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); } } - + tmp *= -1.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + 2 * source_fe.degree, j + (k + 4) * this->degree - k - 2) = tmp; @@ -1227,44 +1223,44 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, } } } - + break; } - + case 3: { for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_face; ++vertex) interpolation_matrix (0, vertex) = 0.25; - + for (unsigned int i = 0; i < this->degree - 1;) { for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_face; ++line) interpolation_matrix (0, i + line * (this->degree - 1) + 4) = -0.5; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (0, i + (j + 4) * this->degree - j) = 1.0; j = j + 2; } - + interpolation_matrix (1, i + 4) = -1.0; interpolation_matrix (2, i + 3 * this->degree + 1) = -1.0; i = i + 2; } - + interpolation_matrix (1, 0) = 0.5; interpolation_matrix (1, 1) = 0.5; interpolation_matrix (2, 2) = 0.5; interpolation_matrix (2, 3) = 0.5; interpolation_matrix (3, 3) = 1.0; - + int factorial_i = 1; int factorial_ij; int factorial_j; int factorial_k; int factorial_kl; int factorial_l; - + for (int i = 2; i <= (int) this->degree; ++i) { double tmp = std::pow (0.5, i + 1); @@ -1273,26 +1269,26 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + 2 * source_fe.degree, i + 2 * this->degree) = tmp; interpolation_matrix (i + 2 * source_fe.degree, i + 3 * this->degree - 1) = tmp; tmp *= -2.0; - + for (unsigned int j = 0; j < this->degree - 1;) { interpolation_matrix (i + 2, j + (i + 2) * this->degree + 2 - i) = tmp; interpolation_matrix (i + 2 * source_fe.degree, i + (j + 4) * this->degree - 2) = tmp; j = j + 2; } - + tmp *= -1.0; interpolation_matrix (i + source_fe.degree + 1, i + this->degree + 1) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, i + 3 * this->degree - 1) = tmp; factorial_k = 1; - + for (int j = 2; j <= (int) this->degree; ++j) { interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (j + 2) * this->degree - j) = std::pow (0.5, i + j); factorial_k *= j; factorial_l = factorial_k; factorial_kl = 1; - + for (int k = j + 1; k <= (int) this->degree; ++k) { factorial_kl *= k - j; @@ -1300,11 +1296,11 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + (j + 2) * source_fe.degree - j, i + (k + 2) * this->degree - k) = std::pow (0.5, i + k) * factorial_l / (factorial_k * factorial_kl); } } - + factorial_i *= i; factorial_j = factorial_i; factorial_ij = 1; - + for (int j = i + 1; j <= (int) this->degree; ++j) { factorial_ij *= j - i; @@ -1318,14 +1314,14 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + source_fe.degree + 1, j + this->degree + 1) = tmp; interpolation_matrix (i + 3 * source_fe.degree - 1, j + 3 * this->degree - 1) = tmp; factorial_k = 1; - + for (int k = 2; k <= (int) this->degree; ++k) { interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (k + 2) * this->degree - k) = tmp * std::pow (0.5, k); factorial_k *= k; factorial_l = factorial_k; factorial_kl = 1; - + for (int l = k + 1; l <= (int) this->degree; ++l) { factorial_kl *= l - k; @@ -1333,15 +1329,15 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, interpolation_matrix (i + (k + 2) * source_fe.degree - k, j + (l + 2) * this->degree - l) = tmp * std::pow (0.5, l) * factorial_l / (factorial_k * factorial_kl); } } - + tmp *= -1.0; - + for (unsigned int k = 0; k < this->degree - 1;) { interpolation_matrix (i + 2, k + (j + 2) * this->degree + 2 - j) = tmp; interpolation_matrix (i + 2 * source_fe.degree, j + (k + 4) * this->degree - 2) = tmp; k = k + 2; - } + } } } } @@ -1350,23 +1346,22 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe, } } -#endif template void FE_Q_Hierarchical::initialize_unit_face_support_points () { const unsigned int codim = dim-1; - + // number of points: (degree+1)^codim unsigned int n = this->degree+1; for (unsigned int i=1; idegree+1; - + this->unit_face_support_points.resize(n); - + Point p; - + unsigned int k=0; for (unsigned int iz=0; iz <= ((codim>2) ? this->degree : 0) ; ++iz) for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy) @@ -1577,9 +1572,9 @@ hierarchic_to_fe_q_hierarchical_numbering (const FiniteElementData &fe) h2l[next_index++] = (2+i)*n2+(2+j)*n+2+k; Assert (next_index == fe.dofs_per_cell, ExcInternalError()); - + break; - } + } default: Assert (false, ExcNotImplemented()); @@ -1599,7 +1594,6 @@ face_fe_q_hierarchical_to_hierarchic_numbering (const unsigned int degree) } -#if (deal_II_dimension == 1) template <> std::vector @@ -1630,7 +1624,6 @@ FE_Q_Hierarchical<1>::has_support_on_face (const unsigned int shape_index, ((shape_index == 1) && (face_index == 1))); } -#endif @@ -1652,10 +1645,10 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, || ((dim==3) && (shape_index>=this->first_hex_index))) return false; - + // let's see whether this is a // vertex - if (shape_index < this->first_line_index) + if (shape_index < this->first_line_index) { // for Q elements, there is // one dof per vertex, so @@ -1677,7 +1670,7 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, = (shape_index - this->first_line_index) / this->dofs_per_line; Assert (line_index < GeometryInfo::lines_per_cell, ExcInternalError()); - + for (unsigned int i=0; i::lines_per_face; ++i) if (GeometryInfo::face_to_cell_lines(face_index,i) == line_index) return true; @@ -1686,18 +1679,18 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, else if (shape_index < this->first_hex_index) // dof is on a quad { - const unsigned int quad_index + const unsigned int quad_index = (shape_index - this->first_quad_index) / this->dofs_per_quad; Assert (static_cast(quad_index) < static_cast(GeometryInfo::quads_per_cell), ExcInternalError()); - + // in 2d, cell bubble are // zero on all faces. but // we have treated this // case above already Assert (dim != 2, ExcInternalError()); - + // in 3d, // quad_index=face_index if (dim == 3) @@ -1723,7 +1716,7 @@ FE_Q_Hierarchical::has_support_on_face (const unsigned int shape_index, template -std::vector +std::vector FE_Q_Hierarchical::get_embedding_dofs (const unsigned int sub_degree) const { Assert ((sub_degree>0) && (sub_degree<=this->degree), @@ -1737,7 +1730,7 @@ FE_Q_Hierarchical::get_embedding_dofs (const unsigned int sub_degree) const return embedding_dofs; } - + if (sub_degree==1) { std::vector embedding_dofs (GeometryInfo::vertices_per_cell); @@ -1757,11 +1750,11 @@ FE_Q_Hierarchical::get_embedding_dofs (const unsigned int sub_degree) const if ((dim==2) || (dim==3)) { - std::vector embedding_dofs ( (dim==2) ? + std::vector embedding_dofs ( (dim==2) ? (sub_degree+1) * (sub_degree+1) : (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); - for (unsigned int i=0; i<( (dim==2) ? + for (unsigned int i=0; i<( (dim==2) ? (sub_degree+1) * (sub_degree+1) : (sub_degree+1) * (sub_degree+1) * (sub_degree+1) ); ++i) { @@ -1769,36 +1762,36 @@ FE_Q_Hierarchical::get_embedding_dofs (const unsigned int sub_degree) const if (i::vertices_per_cell) embedding_dofs[i] = i; // line - else if (i<(GeometryInfo::vertices_per_cell + + else if (i<(GeometryInfo::vertices_per_cell + GeometryInfo::lines_per_cell * (sub_degree-1))) { const unsigned int j = (i - GeometryInfo::vertices_per_cell) % (sub_degree-1); const unsigned int line = (i - GeometryInfo::vertices_per_cell - j) / (sub_degree-1); - embedding_dofs[i] = GeometryInfo::vertices_per_cell + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + line * (this->degree-1) + j; } // quad - else if (i<(GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (sub_degree-1)) + + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1)) + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1)) { const unsigned int j = (i - GeometryInfo::vertices_per_cell - GeometryInfo::lines_per_cell * (sub_degree-1)) % (sub_degree-1); const unsigned int k = ( (i - GeometryInfo::vertices_per_cell - GeometryInfo::lines_per_cell * (sub_degree-1) - j) / (sub_degree-1) ) % (sub_degree-1); - const unsigned int face = (i - GeometryInfo::vertices_per_cell - + const unsigned int face = (i - GeometryInfo::vertices_per_cell - GeometryInfo::lines_per_cell * (sub_degree-1) - k * (sub_degree-1) - j) / ( (sub_degree-1) * (sub_degree-1) ); - embedding_dofs[i] = GeometryInfo::vertices_per_cell + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + GeometryInfo::lines_per_cell * (this->degree-1) + face * (this->degree-1) * (this->degree-1) + k * (this->degree-1) + j; } // hex - else if (i<(GeometryInfo::vertices_per_cell + - GeometryInfo::lines_per_cell * (sub_degree-1)) + + else if (i<(GeometryInfo::vertices_per_cell + + GeometryInfo::lines_per_cell * (sub_degree-1)) + GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) + GeometryInfo::hexes_per_cell * (sub_degree-1) * (sub_degree-1) * (sub_degree-1)) { @@ -1811,8 +1804,8 @@ FE_Q_Hierarchical::get_embedding_dofs (const unsigned int sub_degree) const const unsigned int l = (i - GeometryInfo::vertices_per_cell - GeometryInfo::lines_per_cell * (sub_degree-1) - GeometryInfo::quads_per_cell * (sub_degree-1) * (sub_degree-1) - j - k * (sub_degree-1)) / ( (sub_degree-1) * (sub_degree-1) ); - - embedding_dofs[i] = GeometryInfo::vertices_per_cell + + + embedding_dofs[i] = GeometryInfo::vertices_per_cell + GeometryInfo::lines_per_cell * (this->degree-1) + GeometryInfo::quads_per_cell * (this->degree-1) * (this->degree-1) + l * (this->degree-1) * (this->degree-1) + k * (this->degree-1) + j; @@ -1840,6 +1833,8 @@ FE_Q_Hierarchical::memory_consumption () const -template class FE_Q_Hierarchical; +// explicit instantiations +#include "fe_q_hierarchical.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_q_hierarchical.inst.in b/deal.II/source/fe/fe_q_hierarchical.inst.in new file mode 100644 index 0000000000..0350400b9f --- /dev/null +++ b/deal.II/source/fe/fe_q_hierarchical.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_Q_Hierarchical; + } + diff --git a/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/source/fe/fe_raviart_thomas.cc index d9a6f3dfdf..71162802a8 100644 --- a/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/source/fe/fe_raviart_thomas.cc @@ -49,7 +49,7 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) { Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - + this->mapping_type = mapping_raviart_thomas; // First, initialize the // generalized support points and @@ -60,7 +60,7 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) //matrix, generating the correct //basis functions from the raw //ones. - + // We use an auxiliary matrix in // this function. Therefore, // inverse_node_matrix is still @@ -73,7 +73,7 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int deg) // From now on, the shape functions // will be the correct ones, not // the raw shape functions anymore. - + // Reinit the vectors of // restriction and prolongation // matrices to the right sizes. @@ -118,7 +118,7 @@ FE_RaviartThomas::get_name () const // polynomial degree and is thus one higher // than the argument given to the // constructor - std::ostringstream namebuf; + std::ostringstream namebuf; namebuf << "FE_RaviartThomas<" << dim << ">(" << this->degree-1 << ")"; return namebuf.str(); @@ -138,58 +138,24 @@ FE_RaviartThomas::clone() const //--------------------------------------------------------------------------- -#if deal_II_dimension == 1 - template void FE_RaviartThomas::initialize_support_points (const unsigned int deg) { - return; - - Assert (false, ExcNotImplemented()); - QGauss cell_quadrature(deg+1); const unsigned int n_interior_points = (deg>0) ? cell_quadrature.size() : 0; - - this->generalized_support_points.resize (2 + n_interior_points); - - // Number of the point being entered - unsigned int current = 0; - - - if (deg==0) return; - interior_weights.reinit(TableIndices<3>(2+n_interior_points, 0, dim)); - - for (unsigned int k=0;kgeneralized_support_points[current++] = cell_quadrature.point(k); - - Assert (current == this->generalized_support_points.size(), - ExcInternalError()); -} - -#else - -// Version for 2d and higher. See above for 1d version -template -void -FE_RaviartThomas::initialize_support_points (const unsigned int deg) -{ - QGauss cell_quadrature(deg+1); - const unsigned int n_interior_points - = (deg>0) ? cell_quadrature.size() : 0; - unsigned int n_face_points = (dim>1) ? 1 : 0; // compute (deg+1)^(dim-1) for (unsigned int d=1;dgeneralized_support_points.resize (GeometryInfo::faces_per_cell*n_face_points + n_interior_points); this->generalized_face_support_points.resize (n_face_points); - + // Number of the point being entered unsigned int current = 0; @@ -200,10 +166,10 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) = Polynomials::Legendre::generate_complete_basis(deg); boundary_weights.reinit(n_face_points, legendre.n()); - + // Assert (face_points.size() == this->dofs_per_face, // ExcInternalError()); - + for (unsigned int k=0;kgeneralized_face_support_points[k] = face_points.point(k); @@ -227,9 +193,9 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) this->generalized_support_points[current] = faces.point(current+QProjector::DataSetDescriptor::face(0,true,false,false,n_face_points)); } } - + if (deg==0) return; - + // Create Legendre basis for the // space D_xi Q_k std::vector* > polynomials(dim); @@ -242,9 +208,9 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) polynomials[dd] = new AnisotropicPolynomials(poly); } - + interior_weights.reinit(TableIndices<3>(n_interior_points, polynomials[0]->n(), dim)); - + for (unsigned int k=0;kgeneralized_support_points[current++] = cell_quadrature.point(k); @@ -256,15 +222,12 @@ FE_RaviartThomas::initialize_support_points (const unsigned int deg) for (unsigned int d=0;dgeneralized_support_points.size(), ExcInternalError()); } -#endif - -#if deal_II_dimension == 1 template <> void @@ -277,7 +240,7 @@ FE_RaviartThomas<1>::initialize_restriction() this->restriction[0][i].reinit(0,0); } -#endif + // This function is the same Raviart-Thomas interpolation performed by // interpolate. Still, we cannot use interpolate, since it was written @@ -328,7 +291,7 @@ FE_RaviartThomas::initialize_restriction() const unsigned int child = GeometryInfo::child_cell_on_face( RefinementCase::isotropic_refinement, face, sub); - + // On a certain face, we must // compute the moments of ALL // fine level functions with @@ -360,9 +323,9 @@ FE_RaviartThomas::initialize_restriction() } } } - + if (this->degree == 1) return; - + // Create Legendre basis for the // space D_xi Q_k. Here, we cannot // use the shape functions @@ -376,7 +339,7 @@ FE_RaviartThomas::initialize_restriction() polynomials[dd] = new AnisotropicPolynomials(poly); } - + QGauss q_cell(this->degree); const unsigned int start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; @@ -389,11 +352,11 @@ FE_RaviartThomas::initialize_restriction() for (unsigned int i = 0; i < this->dofs_per_cell; ++i) for (unsigned int d=0;dshape_value_component(i, q_cell.point(k), d); - + for (unsigned int child=0;child::max_children_per_cell;++child) { Quadrature q_sub = QProjector::project_to_child(q_cell, child); - + for (unsigned int k=0;kdofs_per_cell; ++i_child) for (unsigned int d=0;d::initialize_restriction() * polynomials[d]->compute_value(i_weight, q_sub.point(k)); } } - + for (unsigned int d=0;d -std::vector -FE_RaviartThomas<1>::get_dpo_vector (const unsigned int deg) -{ - std::vector dpo(2); - dpo[0] = 1; - dpo[1] = deg; - return dpo; -} - -#endif - template std::vector @@ -440,11 +389,11 @@ FE_RaviartThomas::get_dpo_vector (const unsigned int deg) // and then there are interior dofs const unsigned int interior_dofs = dim*deg*dofs_per_face; - + std::vector dpo(dim+1); dpo[dim-1] = dofs_per_face; dpo[dim] = interior_dofs; - + return dpo; } @@ -484,16 +433,16 @@ FE_RaviartThomas::has_support_on_face ( // these in a table return (face_index != GeometryInfo::opposite_face[shape_index]); } - + default: return true; }; }; - + default: // other rt_order return true; }; - + return true; } @@ -534,10 +483,10 @@ FE_RaviartThomas::interpolate( local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i) * values[face*n_face_points+k](GeometryInfo::unit_normal_direction[face]+offset); } - + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; - + for (unsigned int k=0;k::interpolate( local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i) * values[GeometryInfo::unit_normal_direction[face]][face*n_face_points+k]; } - + const unsigned start_cell_dofs = GeometryInfo::faces_per_cell*this->dofs_per_face; const unsigned start_cell_points = GeometryInfo::faces_per_cell*n_face_points; - + for (unsigned int k=0;k::memory_consumption () const } -template class FE_RaviartThomas; + +// explicit instantiations +#include "fe_raviart_thomas.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_raviart_thomas.inst.in b/deal.II/source/fe/fe_raviart_thomas.inst.in new file mode 100644 index 0000000000..46eca9e9b1 --- /dev/null +++ b/deal.II/source/fe/fe_raviart_thomas.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_RaviartThomas; + } + diff --git a/deal.II/source/fe/fe_raviart_thomas_nodal.cc b/deal.II/source/fe/fe_raviart_thomas_nodal.cc index 10221ca74c..716d9965b0 100644 --- a/deal.II/source/fe/fe_raviart_thomas_nodal.cc +++ b/deal.II/source/fe/fe_raviart_thomas_nodal.cc @@ -42,7 +42,7 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) { Assert (dim >= 2, ExcImpossibleInDim(dim)); const unsigned int n_dofs = this->dofs_per_cell; - + this->mapping_type = mapping_raviart_thomas; // First, initialize the // generalized support points and @@ -53,7 +53,7 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) //matrix, generating the correct //basis functions from the raw //ones. - + // We use an auxiliary matrix in // this function. Therefore, // inverse_node_matrix is still @@ -66,7 +66,7 @@ FE_RaviartThomasNodal::FE_RaviartThomasNodal (const unsigned int deg) // From now on, the shape functions // will be the correct ones, not // the raw shape functions anymore. - + // Reinit the vectors of // prolongation matrices to the // right sizes. There are no @@ -115,7 +115,7 @@ FE_RaviartThomasNodal::get_name () const // polynomial degree and is thus one higher // than the argument given to the // constructor - std::ostringstream namebuf; + std::ostringstream namebuf; namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << this->degree-1 << ")"; return namebuf.str(); @@ -172,7 +172,7 @@ FE_RaviartThomasNodal::initialize_support_points (const unsigned int deg) current = this->dofs_per_face*GeometryInfo::faces_per_cell; } - + if (deg==0) return; // In the interior, we need // anisotropic Gauss quadratures, @@ -190,7 +190,7 @@ FE_RaviartThomasNodal::initialize_support_points (const unsigned int deg) ((d==1) ? low : high), ((d==2) ? low : high)); Assert(dim<=3, ExcNotImplemented()); - + for (unsigned int k=0;ksize();++k) this->generalized_support_points[current++] = quadrature->point(k); delete quadrature; @@ -199,20 +199,6 @@ FE_RaviartThomasNodal::initialize_support_points (const unsigned int deg) } -#if deal_II_dimension == 1 - -template <> -std::vector -FE_RaviartThomasNodal<1>::get_dpo_vector (const unsigned int deg) -{ - std::vector dpo(2); - dpo[0] = 1; - dpo[1] = deg; - return dpo; -} - -#endif - template std::vector @@ -227,18 +213,16 @@ FE_RaviartThomasNodal::get_dpo_vector (const unsigned int deg) // and then there are interior dofs const unsigned int interior_dofs = dim*deg*dofs_per_face; - + std::vector dpo(dim+1); dpo[dim-1] = dofs_per_face; dpo[dim] = interior_dofs; - + return dpo; } -#if deal_II_dimension == 1 - template <> std::vector FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int) @@ -247,7 +231,6 @@ FE_RaviartThomasNodal<1>::get_ria_vector (const unsigned int) return std::vector(); } -#endif template @@ -287,14 +270,14 @@ FE_RaviartThomasNodal::has_support_on_face ( // on the faces and each face has // degree degrees. const unsigned int support_face = shape_index / this->degree; - + // The only thing we know for sure // is that shape functions with // support on one face are zero on // the opposite face. if (support_face < GeometryInfo::faces_per_cell) return (face_index != GeometryInfo::opposite_face[support_face]); - + // In all other cases, return true, // which is safe return true; @@ -429,7 +412,7 @@ FE_RaviartThomasNodal::hp_vertex_dof_identities ( else { Assert(false, ExcNotImplemented()); - return std::vector > (); + return std::vector > (); } } @@ -450,7 +433,7 @@ hp_line_dof_identities (const FiniteElement &fe_other) const // only lines in 2d if (dim != 2) return std::vector >(); - + // dofs are located along lines, so two // dofs are identical only if in the // following two cases (remember that @@ -472,7 +455,7 @@ hp_line_dof_identities (const FiniteElement &fe_other) const // the constructor const unsigned int p = this->degree-1; const unsigned int q = fe_q_other->degree-1; - + std::vector > identities; if (p==q) @@ -482,7 +465,7 @@ hp_line_dof_identities (const FiniteElement &fe_other) const else if (p%2==0 && q%2==0) identities.push_back(std::make_pair(p/2,q/2)); - + return identities; } else @@ -513,7 +496,7 @@ FE_RaviartThomasNodal::hp_quad_dof_identities ( // case above const unsigned int p = this->dofs_per_quad; const unsigned int q = fe_q_other->dofs_per_quad; - + std::vector > identities; if (p==q) @@ -549,13 +532,12 @@ FE_RaviartThomasNodal::compare_for_face_domination ( else return FiniteElementDomination::other_element_dominates; } - + Assert (false, ExcNotImplemented()); return FiniteElementDomination::neither_element_dominates; } -#if deal_II_dimension == 1 template <> void @@ -577,11 +559,8 @@ FE_RaviartThomasNodal<1>::get_subface_interpolation_matrix ( Assert (false, ExcImpossibleInDim(1)); } -#endif -#if deal_II_dimension > 1 - template void FE_RaviartThomasNodal::get_face_interpolation_matrix ( @@ -596,7 +575,7 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix ( (dynamic_cast*>(&x_source_fe) != 0), typename FiniteElement:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.n() == this->dofs_per_face, ExcDimensionMismatch (interpolation_matrix.n(), this->dofs_per_face)); @@ -623,7 +602,7 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix ( Assert (this->dofs_per_face <= source_fe.dofs_per_face, typename FiniteElement:: ExcInterpolationNotImplemented ()); - + // generate a quadrature // with the generalized support points. // This is later based as a @@ -631,26 +610,26 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix ( // which returns the support // points on the face. Quadrature quad_face_support (source_fe.get_generalized_face_support_points ()); - + // Rule of thumb for FP accuracy, // that can be expected for a // given polynomial degree. // This value is used to cut // off values close to zero. double eps = 2e-13*this->degree*(dim-1); - + // compute the interpolation // matrix by simply taking the // value at the support points. const Quadrature face_projection = QProjector::project_to_face (quad_face_support, 0); - + for (unsigned int i=0; i &p = face_projection.point (i); for (unsigned int j=0; jdofs_per_face; ++j) - { + { double matrix_entry = this->shape_value_component (this->face_to_cell_index(j, 0), p, 0); @@ -667,7 +646,7 @@ FE_RaviartThomasNodal::get_face_interpolation_matrix ( matrix_entry = 0.0; interpolation_matrix(i,j) = matrix_entry; - } + } } // make sure that the row sum of @@ -703,7 +682,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix ( (dynamic_cast*>(&x_source_fe) != 0), typename FiniteElement:: ExcInterpolationNotImplemented()); - + Assert (interpolation_matrix.n() == this->dofs_per_face, ExcDimensionMismatch (interpolation_matrix.n(), this->dofs_per_face)); @@ -730,7 +709,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix ( Assert (this->dofs_per_face <= source_fe.dofs_per_face, typename FiniteElement:: ExcInterpolationNotImplemented ()); - + // generate a quadrature // with the generalized support points. // This is later based as a @@ -738,14 +717,14 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix ( // which returns the support // points on the face. Quadrature quad_face_support (source_fe.get_generalized_face_support_points ()); - + // Rule of thumb for FP accuracy, // that can be expected for a // given polynomial degree. // This value is used to cut // off values close to zero. double eps = 2e-13*this->degree*(dim-1); - + // compute the interpolation // matrix by simply taking the // value at the support points. @@ -758,7 +737,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix ( const Point &p = subface_projection.point (i); for (unsigned int j=0; jdofs_per_face; ++j) - { + { double matrix_entry = this->shape_value_component (this->face_to_cell_index(j, 0), p, 0); @@ -774,7 +753,7 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix ( matrix_entry = 0.0; interpolation_matrix(i,j) = matrix_entry; - } + } } // make sure that the row sum of @@ -794,10 +773,10 @@ FE_RaviartThomasNodal::get_subface_interpolation_matrix ( } } -#endif +// explicit instantiations +#include "fe_raviart_thomas_nodal.inst" -template class FE_RaviartThomasNodal; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_raviart_thomas_nodal.inst.in b/deal.II/source/fe/fe_raviart_thomas_nodal.inst.in new file mode 100644 index 0000000000..951a219fe4 --- /dev/null +++ b/deal.II/source/fe/fe_raviart_thomas_nodal.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FE_RaviartThomasNodal; + } + diff --git a/deal.II/source/fe/fe_system.cc b/deal.II/source/fe/fe_system.cc index 6090ca2a82..098a75e02f 100644 --- a/deal.II/source/fe/fe_system.cc +++ b/deal.II/source/fe/fe_system.cc @@ -1950,33 +1950,16 @@ initialize_unit_support_points () } -#if deal_II_dimension == 1 - -template <> -void -FESystem<1>:: -initialize_unit_face_support_points () -{ - // no faces no work -} - - -template <> -void -FESystem<1,2>:: -initialize_unit_face_support_points () -{ - // no faces no work -} - -#endif - template void FESystem:: initialize_unit_face_support_points () { + // Nothing to do in 1D + if (dim == 1) + return; + // if one of the base elements has // no support points, then it makes // no sense to define support @@ -2018,21 +2001,15 @@ template void FESystem::initialize_quad_dof_index_permutation () { - // general template for 1D and 2D, do nothing -} - - + // nothing to do in other dimensions than 3 + if (dim < 3) + return; -#if deal_II_dimension == 3 - -template <> -void -FESystem<3>::initialize_quad_dof_index_permutation () -{ // the array into which we want to write // should have the correct size already. - Assert (adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad, - ExcInternalError()); + Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()== + 8*this->dofs_per_quad, ExcInternalError()); + // to obtain the shifts for this composed // element, copy the shift information of // the base elements @@ -2045,7 +2022,8 @@ FESystem<3>::initialize_quad_dof_index_permutation () { for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation_table(index+i,j)= + temp(i,j); index += temp.size(0); } } @@ -2054,8 +2032,8 @@ FESystem<3>::initialize_quad_dof_index_permutation () // aditionally compose the permutation // information for lines - Assert (adjust_line_dof_index_for_line_orientation_table.size()==this->dofs_per_line, - ExcInternalError()); + Assert (this->adjust_line_dof_index_for_line_orientation_table.size()== + this->dofs_per_line, ExcInternalError()); index = 0; for (unsigned int b=0; b::initialize_quad_dof_index_permutation () for (unsigned int c=0; cadjust_line_dof_index_for_line_orientation_table.begin() + +index); index += temp2.size(); } } @@ -2072,8 +2051,6 @@ FESystem<3>::initialize_quad_dof_index_permutation () ExcInternalError()); } -#endif - template @@ -3407,9 +3384,6 @@ FESystem::memory_consumption () const // explicit instantiations -template class FESystem; -#if deal_II_dimension != 3 -template class FESystem; -#endif +#include "fe_system.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_system.inst.in b/deal.II/source/fe/fe_system.inst.in new file mode 100644 index 0000000000..97153d8622 --- /dev/null +++ b/deal.II/source/fe/fe_system.inst.in @@ -0,0 +1,22 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class FESystem; +#if deal_II_dimension != 3 + template class FESystem; +#endif + } + diff --git a/deal.II/source/fe/fe_tools.all_dimensions.cc b/deal.II/source/fe/fe_tools.all_dimensions.cc deleted file mode 100644 index 4c40c955fb..0000000000 --- a/deal.II/source/fe/fe_tools.all_dimensions.cc +++ /dev/null @@ -1,302 +0,0 @@ -//--------------------------------------------------------------------------- -// $Id$ -// Version: $Name$ -// -// Copyright (C) 2005, 2006, 2008, 2010 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. -// -//--------------------------------------------------------------------------- - - -#include -#include -#include - -DEAL_II_NAMESPACE_OPEN - -namespace FETools -{ - - template - void - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, - std::vector &h2l) - { - Assert (h2l.size() == fe.dofs_per_cell, - ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); - h2l = hierarchic_to_lexicographic_numbering (fe); - } - - - - template - std::vector - hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) - { - Assert (fe.n_components() == 1, ExcInvalidFE()); - - std::vector h2l (fe.dofs_per_cell); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - // polynomial degree - const unsigned int degree = fe.dofs_per_line+1; - // number of grid points in each - // direction - const unsigned int n = degree+1; - - // the following lines of code are - // somewhat odd, due to the way the - // hierarchic numbering is - // organized. if someone would - // really want to understand these - // lines, you better draw some - // pictures where you indicate the - // indices and orders of vertices, - // lines, etc, along with the - // numbers of the degrees of - // freedom in hierarchical and - // lexicographical order - switch (dim) - { - case 1: - { - h2l[0] = 0; - h2l[1] = dofs_per_cell-1; - for (unsigned int i=2; i - void - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, - std::vector &l2h) - { - l2h = lexicographic_to_hierarchic_numbering (fe); - } - - - - template - std::vector - lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) - { - return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); - } - -} - - -// explicit instantiations -namespace FETools -{ - template - void - hierarchic_to_lexicographic_numbering<1> - (const FiniteElementData<1> &fe, - std::vector &h2l); - template - void - hierarchic_to_lexicographic_numbering<2> - (const FiniteElementData<2> &fe, - std::vector &h2l); - template - void - hierarchic_to_lexicographic_numbering<3> - (const FiniteElementData<3> &fe, - std::vector &h2l); - - - template - void - lexicographic_to_hierarchic_numbering<1> - (const FiniteElementData<1> &fe, - std::vector &l2h); - template - void - lexicographic_to_hierarchic_numbering<2> - (const FiniteElementData<2> &fe, - std::vector &l2h); - template - void - lexicographic_to_hierarchic_numbering<3> - (const FiniteElementData<3> &fe, - std::vector &l2h); - - - - template - std::vector - hierarchic_to_lexicographic_numbering<1> - (const FiniteElementData<1> &fe); - template - std::vector - hierarchic_to_lexicographic_numbering<2> - (const FiniteElementData<2> &fe); - template - std::vector - hierarchic_to_lexicographic_numbering<3> - (const FiniteElementData<3> &fe); - - - template - std::vector - lexicographic_to_hierarchic_numbering<1> - (const FiniteElementData<1> &fe); - template - std::vector - lexicographic_to_hierarchic_numbering<2> - (const FiniteElementData<2> &fe); - template - std::vector - lexicographic_to_hierarchic_numbering<3> - (const FiniteElementData<3> &fe); -} - - -DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_tools.cc b/deal.II/source/fe/fe_tools.cc index ec1298b4fc..ff5dfe551a 100644 --- a/deal.II/source/fe/fe_tools.cc +++ b/deal.II/source/fe/fe_tools.cc @@ -53,51 +53,52 @@ DEAL_II_NAMESPACE_OPEN namespace { - // a shared pointer to factory + // use a shared pointer to factory // objects, to ensure that they get // deleted at the end of the // program run and don't end up as // apparent memory leaks to - // programs like valgrind - typedef - std_cxx1x::shared_ptr > - FEFactoryPointer; - + // programs like valgrind. // a function that returns the // default set of finite element // names and factory objects for // them. used to initialize // fe_name_map below + template #ifdef DEAL_II_ANON_NAMESPACE_BUG static #endif - std::map + std::map > > get_default_fe_names () { - std::map default_map; + std::map > > default_map; + + typedef std_cxx1x::shared_ptr > + FEFactoryPointer; default_map["FE_Q_Hierarchical"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_ABF"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_RaviartThomas"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_RaviartThomasNodal"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_Nedelec"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_DGPNonparametric"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_DGP"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_DGPMonomial"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_DGQ"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_DGQArbitraryNodes"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); default_map["FE_Q"] - = FEFactoryPointer(new FETools::FEFactory >); + = FEFactoryPointer(new FETools::FEFactory >); return default_map; } @@ -134,17 +135,29 @@ namespace // static so that we can link // several dimensions together. // - // it is initialized at program - // start time using the function - // above. because at this time - // there are no threads running, - // there are no thread-safety - // issues here + // it is initialized at program start time + // using the function above. because at + // this time there are no threads running, + // there are no thread-safety issues + // here. since this is compiled for all + // dimensions at once, need to create + // objects for each dimension and then + // separate between them further down + static + std::map > > + fe_name_map_1d + = get_default_fe_names<1> (); static std::map > > - fe_name_map - = get_default_fe_names (); + std_cxx1x::shared_ptr > > + fe_name_map_2d + = get_default_fe_names<2> (); + static + std::map > > + fe_name_map_3d + = get_default_fe_names<3> (); } @@ -577,7 +590,8 @@ namespace FETools } } -#if deal_II_dimension == 1 + + template<> void compute_embedding_matrices(const FiniteElement<1,2> &, @@ -587,7 +601,8 @@ namespace FETools Assert(false, ExcNotImplemented()); } -#elif deal_II_dimension == 2 + + template<> void compute_embedding_matrices(const FiniteElement<2,3>&, @@ -597,7 +612,6 @@ namespace FETools Assert(false, ExcNotImplemented()); } -#endif namespace { @@ -763,7 +777,7 @@ namespace FETools } -// This function is tested by tests/fe/internals, since it produces the matrices printed there +// This function is tested by tests/fe/internals, since it produces the matrices printed there template void compute_face_embedding_matrices(const FiniteElement& fe, @@ -773,21 +787,21 @@ namespace FETools { Assert(face_coarse==0, ExcNotImplemented()); Assert(face_fine==0, ExcNotImplemented()); - + const unsigned int nc = GeometryInfo::max_children_per_face; const unsigned int n = fe.dofs_per_face; const unsigned int nd = fe.n_components(); const unsigned int degree = fe.degree; - + const bool normal = fe.conforms(FiniteElementData::Hdiv); const bool tangential = fe.conforms(FiniteElementData::Hcurl); - + for (unsigned int i=0;i H(A); - + Vector v_coarse(nq*nd); Vector v_fine(n); @@ -959,7 +973,8 @@ namespace FETools } } -#if deal_II_dimension == 1 + + template <> void compute_projection_matrices(const FiniteElement<1,2>&, @@ -968,7 +983,8 @@ namespace FETools Assert(false, ExcNotImplemented()); } -#elif deal_II_dimension == 2 + + template <> void compute_projection_matrices(const FiniteElement<2,3>&, @@ -977,7 +993,6 @@ namespace FETools Assert(false, ExcNotImplemented()); } -#endif template @@ -1635,10 +1650,69 @@ for (; cell!=endc; ++cell) } - template + template <> + void + add_fe_name(const std::string& parameter_name, + const FEFactoryBase<1,1>* factory) + { + // Erase everything after the + // actual class name + std::string name = parameter_name; + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + if (name_end < name.size()) + name.erase(name_end); + // first make sure that no other + // thread intercepts the + // operation of this function; + // for this, acquire the lock + // until we quit this function + Threads::ThreadMutex::ScopedLock lock(fe_name_map_lock); + + Assert(fe_name_map_1d.find(name) == fe_name_map_1d.end(), + ExcMessage("Cannot change existing element in finite element name list")); + + // Insert the normalized name into + // the map + fe_name_map_1d[name] = + std_cxx1x::shared_ptr > (factory); + } + + + + template <> + void + add_fe_name(const std::string& parameter_name, + const FEFactoryBase<2,2>* factory) + { + // Erase everything after the + // actual class name + std::string name = parameter_name; + unsigned int name_end = + name.find_first_not_of(std::string("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_")); + if (name_end < name.size()) + name.erase(name_end); + // first make sure that no other + // thread intercepts the + // operation of this function; + // for this, acquire the lock + // until we quit this function + Threads::ThreadMutex::ScopedLock lock(fe_name_map_lock); + + Assert(fe_name_map_2d.find(name) == fe_name_map_2d.end(), + ExcMessage("Cannot change existing element in finite element name list")); + + // Insert the normalized name into + // the map + fe_name_map_2d[name] = + std_cxx1x::shared_ptr > (factory); + } + + + template <> void add_fe_name(const std::string& parameter_name, - const FEFactoryBase* factory) + const FEFactoryBase<3,3>* factory) { // Erase everything after the // actual class name @@ -1654,12 +1728,13 @@ for (; cell!=endc; ++cell) // until we quit this function Threads::ThreadMutex::ScopedLock lock(fe_name_map_lock); - Assert(fe_name_map.find(name) == fe_name_map.end(), + Assert(fe_name_map_3d.find(name) == fe_name_map_3d.end(), ExcMessage("Cannot change existing element in finite element name list")); // Insert the normalized name into // the map - fe_name_map[name] = FEFactoryPointer(factory); + fe_name_map_3d[name] = + std_cxx1x::shared_ptr > (factory); } @@ -1667,9 +1742,16 @@ for (; cell!=endc; ++cell) { namespace { + // TODO: this encapsulates the call to the + // dimension-dependent fe_name_map so that we + // have a unique interface. could be done + // smarter? template FiniteElement* - get_fe_from_name (std::string &name) + get_fe_from_name_ext (std::string &name, + const std::map > > + &fe_name_map) { // Extract the name of the // finite element class, which only @@ -1713,7 +1795,8 @@ for (; cell!=endc; ++cell) // Now, the name of the // first base element is // first... Let's get it - base_fes.push_back (get_fe_from_name (name)); + base_fes.push_back (get_fe_from_name_ext (name, + fe_name_map)); // next check whether // FESystem placed a // multiplicity after @@ -1892,6 +1975,35 @@ for (; cell!=endc; ++cell) // here after throwing return 0; } + + + + // need to work around problem with different + // name map names for different dimensions + // TODO: should be possible to do nicer! + template + FiniteElement* get_fe_from_name (std::string &name); + + template <> + FiniteElement<1,1>* + get_fe_from_name<1> (std::string &name) + { + return get_fe_from_name_ext<1,1> (name, fe_name_map_1d); + } + + template <> + FiniteElement<2,2>* + get_fe_from_name<2> (std::string &name) + { + return get_fe_from_name_ext<2,2> (name, fe_name_map_2d); + } + + template <> + FiniteElement<3,3>* + get_fe_from_name<3> (std::string &name) + { + return get_fe_from_name_ext<3,3> (name, fe_name_map_3d); + } } } @@ -2220,495 +2332,222 @@ for (; cell!=endc; ++cell) Assert (X.n() == rhs_quadrature.size(), ExcInternalError()); } -} + template + void + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe, + std::vector &h2l) + { + Assert (h2l.size() == fe.dofs_per_cell, + ExcDimensionMismatch (h2l.size(), fe.dofs_per_cell)); + h2l = hierarchic_to_lexicographic_numbering (fe); + } -/*-------------- Explicit Instantiations -------------------------------*/ - -namespace FETools -{ - template class FEFactoryBase; - - template - void compute_node_matrix( - FullMatrix&, - const FiniteElement&); - template - void compute_component_wise( - const FiniteElement& element, - std::vector&, std::vector >&); - template - void compute_block_renumbering ( - const FiniteElement& element, - std::vector&, std::vector&_indices, bool); - template - void get_interpolation_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - template - void get_back_interpolation_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - template - void get_interpolation_difference_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - template - void get_interpolation_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - template - void get_back_interpolation_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - template - void get_interpolation_difference_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - - template - void get_projection_matrix - (const FiniteElement &, - const FiniteElement &, - FullMatrix &); - - template - void compute_embedding_matrices - (const FiniteElement &, std::vector > >&,bool); - - template - void compute_face_embedding_matrices - (const FiniteElement &, FullMatrix (&matrices)[GeometryInfo::max_children_per_face], - unsigned int, unsigned int); - - template - void compute_projection_matrices - (const FiniteElement &, std::vector > >&, bool); - - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - -#if deal_II_dimension != 3 - template - void compute_block_renumbering ( - const FiniteElement& element, - std::vector&, std::vector&_indices, bool); - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); -#endif - template - void back_interpolate - (const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); - template - void back_interpolate - (const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - template - void interpolation_difference - (const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); - template - void interpolation_difference - (const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - template - void project_dg - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void extrapolate - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void extrapolate - (const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - - - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - -#if deal_II_dimension != 3 - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void interpolate - (const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); -#endif - - template - void back_interpolate - (const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); - template - void back_interpolate - (const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - template - void interpolation_difference - (const DoFHandler &, const Vector &, - const FiniteElement &, Vector &); - template - void interpolation_difference - (const DoFHandler &, const ConstraintMatrix &, - const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - template - void project_dg - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void extrapolate - (const DoFHandler &, const Vector &, - const DoFHandler &, Vector &); - template - void extrapolate - (const DoFHandler &, const Vector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - - -#ifdef DEAL_II_USE_TRILINOS - - template - void interpolate - (const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, TrilinosWrappers::Vector &); - template - void interpolate - (const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); - template - void back_interpolate - (const DoFHandler &, const TrilinosWrappers::Vector &, - const FiniteElement &, TrilinosWrappers::Vector &); - template - void back_interpolate - (const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); - template - void interpolation_difference - (const DoFHandler &, const TrilinosWrappers::Vector &, - const FiniteElement &, TrilinosWrappers::Vector &); - template - void interpolation_difference - (const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); - template - void project_dg - (const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, TrilinosWrappers::Vector &); - template - void extrapolate - (const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, TrilinosWrappers::Vector &); - template - void extrapolate - (const DoFHandler &, const TrilinosWrappers::Vector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::Vector &); + template + std::vector + hierarchic_to_lexicographic_numbering (const FiniteElementData &fe) + { + Assert (fe.n_components() == 1, ExcInvalidFE()); + + std::vector h2l (fe.dofs_per_cell); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + // polynomial degree + const unsigned int degree = fe.dofs_per_line+1; + // number of grid points in each + // direction + const unsigned int n = degree+1; + + // the following lines of code are + // somewhat odd, due to the way the + // hierarchic numbering is + // organized. if someone would + // really want to understand these + // lines, you better draw some + // pictures where you indicate the + // indices and orders of vertices, + // lines, etc, along with the + // numbers of the degrees of + // freedom in hierarchical and + // lexicographical order + switch (dim) + { + case 1: + { + h2l[0] = 0; + h2l[1] = dofs_per_cell-1; + for (unsigned int i=2; i - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - -#if deal_II_dimension != 3 - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -#endif + case 3: + { + unsigned int next_index = 0; + // first the eight vertices + h2l[next_index++] = 0; // 0 + h2l[next_index++] = ( 1)*degree; // 1 + h2l[next_index++] = ( n )*degree; // 2 + h2l[next_index++] = ( n+1)*degree; // 3 + h2l[next_index++] = (n*n )*degree; // 4 + h2l[next_index++] = (n*n +1)*degree; // 5 + h2l[next_index++] = (n*n+n )*degree; // 6 + h2l[next_index++] = (n*n+n+1)*degree; // 7 + + // line 0 + for (unsigned int i=0; i - (const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); - template - void back_interpolate - (const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - template - void interpolation_difference - (const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); - template - void interpolation_difference - (const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - template - void project_dg - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, Vector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - - - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - -#if deal_II_dimension != 3 - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void interpolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); -#endif + default: + Assert (false, ExcNotImplemented()); + } - template - void back_interpolate - (const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); - template - void back_interpolate - (const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - template - void interpolation_difference - (const DoFHandler &, const BlockVector &, - const FiniteElement &, BlockVector &); - template - void interpolation_difference - (const DoFHandler &, const ConstraintMatrix &, - const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - template - void project_dg - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, BlockVector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - BlockVector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, Vector &); - template - void extrapolate - (const DoFHandler &, const BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - Vector &); - - -#ifdef DEAL_II_USE_TRILINOS - - template - void interpolate - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, TrilinosWrappers::BlockVector &); - template - void interpolate - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); - template - void back_interpolate - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const FiniteElement &, TrilinosWrappers::BlockVector &); - template - void back_interpolate - (const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); - template - void interpolation_difference - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const FiniteElement &, TrilinosWrappers::BlockVector &); - template - void interpolation_difference - (const DoFHandler &, const ConstraintMatrix &, - const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); - template - void project_dg - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, TrilinosWrappers::BlockVector &); - template - void extrapolate - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, TrilinosWrappers::BlockVector &); - template - void extrapolate - (const DoFHandler &, const TrilinosWrappers::BlockVector &, - const DoFHandler &, const ConstraintMatrix &, - TrilinosWrappers::BlockVector &); + return h2l; + } -#endif - template - void interpolate - (const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, Vector &); - template - void interpolate - (const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, const ConstraintMatrix &, - Vector &); - template - void interpolate - (const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, Vector &); - template - void interpolate - (const hp::DoFHandler &, const Vector &, - const hp::DoFHandler &, const ConstraintMatrix &, - Vector &); - - - template FiniteElement * - get_fe_from_name (const std::string &); - - - template - void add_fe_name( - const std::string& name, - const FEFactoryBase* factory); - - template + template void - compute_projection_from_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - FullMatrix &X); + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe, + std::vector &l2h) + { + l2h = lexicographic_to_hierarchic_numbering (fe); + } - template - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector< Tensor<1, deal_II_dimension > > &vector_of_tensors_at_qp, - std::vector< Tensor<1, deal_II_dimension > > &vector_of_tensors_at_nodes); - template - void - compute_projection_from_quadrature_points( - const FullMatrix &projection_matrix, - const std::vector > &vector_of_tensors_at_qp, - std::vector > &vector_of_tensors_at_nodes); + template + std::vector + lexicographic_to_hierarchic_numbering (const FiniteElementData &fe) + { + return Utilities::invert_permutation(hierarchic_to_lexicographic_numbering (fe)); + } - template - void - compute_interpolation_to_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &quadrature, - FullMatrix &I_q); +} -#if deal_II_dimension != 1 - template - void - compute_projection_from_face_quadrature_points_matrix (const FiniteElement &fe, - const Quadrature &lhs_quadrature, - const Quadrature &rhs_quadrature, - const DoFHandler::active_cell_iterator & cell, - unsigned int face, - FullMatrix &X); -#endif -} +/*-------------- Explicit Instantiations -------------------------------*/ +#include "fe_tools.inst" /*---------------------------- fe_tools.cc ---------------------------*/ diff --git a/deal.II/source/fe/fe_tools.inst.in b/deal.II/source/fe/fe_tools.inst.in new file mode 100644 index 0000000000..16f935646f --- /dev/null +++ b/deal.II/source/fe/fe_tools.inst.in @@ -0,0 +1,522 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + + namespace FETools + \{ + template class FEFactoryBase; + + template + void compute_node_matrix( + FullMatrix&, + const FiniteElement&); + + template + void compute_component_wise( + const FiniteElement& , + std::vector&, std::vector >&); + template + void compute_block_renumbering ( + const FiniteElement& , + std::vector&, std::vector&, bool); + template + void get_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_back_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_interpolation_difference_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_back_interpolation_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + template + void get_interpolation_difference_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + + template + void get_projection_matrix + (const FiniteElement &, + const FiniteElement &, + FullMatrix &); + + template + void compute_embedding_matrices + (const FiniteElement &, std::vector > >&,bool); + + template + void compute_face_embedding_matrices + (const FiniteElement &, FullMatrix (&)[GeometryInfo::max_children_per_face], + unsigned int, unsigned int); + + template + void compute_projection_matrices + (const FiniteElement &, std::vector > >&, bool); + + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + +#if deal_II_dimension != 3 + template + void compute_block_renumbering ( + const FiniteElement& , + std::vector&, std::vector&, bool); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); +#endif + + template + void back_interpolate + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void interpolation_difference + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void project_dg + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + + + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + +#if deal_II_dimension != 3 + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void interpolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); +#endif + + template + void back_interpolate + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void interpolation_difference + (const DoFHandler &, const Vector &, + const FiniteElement &, Vector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void project_dg + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const Vector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + + +#ifdef DEAL_II_USE_TRILINOS + + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, TrilinosWrappers::Vector &); + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + template + void back_interpolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const FiniteElement &, TrilinosWrappers::Vector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + template + void interpolation_difference + (const DoFHandler &, const TrilinosWrappers::Vector &, + const FiniteElement &, TrilinosWrappers::Vector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + template + void project_dg + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, TrilinosWrappers::Vector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, TrilinosWrappers::Vector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::Vector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::Vector &); + +#endif + + + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + +#if deal_II_dimension != 3 + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); +#endif + + template + void back_interpolate + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void project_dg + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + + + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + +#if deal_II_dimension != 3 + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void interpolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); +#endif + + template + void back_interpolate + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const BlockVector &, + const FiniteElement &, BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void project_dg + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + BlockVector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, Vector &); + template + void extrapolate + (const DoFHandler &, const BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + Vector &); + + +#ifdef DEAL_II_USE_TRILINOS + + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, TrilinosWrappers::BlockVector &); + template + void interpolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + template + void back_interpolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const FiniteElement &, TrilinosWrappers::BlockVector &); + template + void back_interpolate + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const FiniteElement &, TrilinosWrappers::BlockVector &); + template + void interpolation_difference + (const DoFHandler &, const ConstraintMatrix &, + const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + template + void project_dg + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, TrilinosWrappers::BlockVector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, TrilinosWrappers::BlockVector &); + template + void extrapolate + (const DoFHandler &, const TrilinosWrappers::BlockVector &, + const DoFHandler &, const ConstraintMatrix &, + TrilinosWrappers::BlockVector &); + +#endif + + + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, Vector &); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, const ConstraintMatrix &, + Vector &); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, Vector &); + template + void interpolate + (const hp::DoFHandler &, const Vector &, + const hp::DoFHandler &, const ConstraintMatrix &, + Vector &); + + + template FiniteElement * + get_fe_from_name (const std::string &); + + + template + void add_fe_name( + const std::string& , + const FEFactoryBase* ); + + template + void + compute_projection_from_quadrature_points_matrix (const FiniteElement &, + const Quadrature &, + const Quadrature &, + FullMatrix &); + + template + void + compute_projection_from_quadrature_points( + const FullMatrix &, + const std::vector< Tensor<1, deal_II_dimension > > &, + std::vector< Tensor<1, deal_II_dimension > > &); + + template + void + compute_projection_from_quadrature_points( + const FullMatrix &, + const std::vector > &, + std::vector > &); + + + template + void + compute_interpolation_to_quadrature_points_matrix (const FiniteElement &, + const Quadrature &, + FullMatrix &); + +#if deal_II_dimension != 1 + template + void + compute_projection_from_face_quadrature_points_matrix (const FiniteElement &, + const Quadrature &, + const Quadrature &, + const DoFHandler::active_cell_iterator & , + unsigned int, + FullMatrix &); +#endif + + template + void + hierarchic_to_lexicographic_numbering + (const FiniteElementData &, + std::vector &); + + template + void + lexicographic_to_hierarchic_numbering + (const FiniteElementData &, + std::vector &); + + template + std::vector + hierarchic_to_lexicographic_numbering + (const FiniteElementData &); + + template + std::vector + lexicographic_to_hierarchic_numbering + (const FiniteElementData &); + \} + } + diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc index 8481733c2f..6727e94170 100644 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@ -3412,32 +3412,16 @@ FEValues::reinit (const typename hp::DoFHandler::cel } -#if deal_II_dimension == 1 - -template <> -void -FEValues<1,2>::reinit (const MGDoFHandler<1,2>::cell_iterator &) -{ - Assert(false, ExcNotImplemented()); -} - -#endif - -#if deal_II_dimension == 2 - -template <> -void -FEValues<2,3>::reinit (const MGDoFHandler<2,3>::cell_iterator &) -{ - Assert(false, ExcNotImplemented()); -} - -#endif - template void FEValues::reinit (const typename MGDoFHandler::cell_iterator &cell) { + if (spacedim > dim) + { + Assert(false, ExcNotImplemented()); + return; + } + // assert that the finite elements // passed to the constructor and // used by the DoFHandler used by @@ -3704,32 +3688,16 @@ void FEFaceValues::reinit (const typename hp::DoFHandler -void FEFaceValues<1,2>::reinit (const MGDoFHandler<1,2>::cell_iterator &, - const unsigned int) -{ - Assert(false,ExcNotImplemented()); -} - -#endif - -#if deal_II_dimension == 2 - -template <> -void FEFaceValues<2,3>::reinit (const MGDoFHandler<2,3>::cell_iterator &, - const unsigned int) -{ - Assert(false,ExcNotImplemented()); -} - -#endif template void FEFaceValues::reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no) { + if (spacedim > dim) + { + Assert(false,ExcNotImplemented()); + return; + } // assert that the finite elements // passed to the constructor and // used by the DoFHandler used by @@ -4137,56 +4105,6 @@ void FESubfaceValues::do_reinit (const unsigned int face_no, /*------------------------------- Explicit Instantiations -------------*/ - -template class FEValuesData; -template class FEValuesBase; -template class FEValues; -template class FEValuesBase:: -CellIterator::cell_iterator>; -template class FEValuesBase:: -CellIterator::cell_iterator>; - -#if deal_II_dimension >= 2 -template class FEFaceValuesBase; -template class FEFaceValues; -template class FESubfaceValues; -#endif - -#if deal_II_dimension != 3 -template class FEValuesData; -template class FEValuesBase; -template class FEValues; -template class FEValuesBase:: -CellIterator::cell_iterator>; -template class FEValuesBase:: - CellIterator::cell_iterator>; - -#if deal_II_dimension == 2 -template class FEFaceValuesBase; -template class FEFaceValues; -template class FESubfaceValues; -#endif -#endif - - -namespace FEValuesViews -{ - template class Scalar; - template class Vector; - template class SymmetricTensor<2, deal_II_dimension, deal_II_dimension>; - -#if deal_II_dimension != 3 - template class Scalar; - template class Vector; - template class SymmetricTensor<2, deal_II_dimension, deal_II_dimension+1>; -#endif -} - - -//--------------------------------------------------------------------------- -// Instantiations are in a different file using the macro IN for the vector type. -// This way, we avoid code reduplication - #include "fe_values.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/fe_values.inst.in b/deal.II/source/fe/fe_values.inst.in index ca6dfc5e1a..502cd78f5b 100644 --- a/deal.II/source/fe/fe_values.inst.in +++ b/deal.II/source/fe/fe_values.inst.in @@ -14,7 +14,55 @@ // Instantiations of functions in FEValuesViews (explicitly write // dealii:: namespace in order to not confuse the compiler with // FEValuesViews::Vector). -for (VEC : SERIAL_VECTORS) +for (deal_II_dimension : DIMENSIONS) +{ +template class FEValuesData; +template class FEValuesBase; +template class FEValues; +template class FEValuesBase:: +CellIterator::cell_iterator>; +template class FEValuesBase:: +CellIterator::cell_iterator>; + +#if deal_II_dimension >= 2 +template class FEFaceValuesBase; +template class FEFaceValues; +template class FESubfaceValues; +#endif + +#if deal_II_dimension != 3 +template class FEValuesData; +template class FEValuesBase; +template class FEValues; +template class FEValuesBase:: +CellIterator::cell_iterator>; +template class FEValuesBase:: + CellIterator::cell_iterator>; + +#if deal_II_dimension == 2 +template class FEFaceValuesBase; +template class FEFaceValues; +template class FESubfaceValues; +#endif +#endif + + +namespace FEValuesViews +\{ + template class Scalar; + template class Vector; + template class SymmetricTensor<2, deal_II_dimension, deal_II_dimension>; + +#if deal_II_dimension != 3 + template class Scalar; + template class Vector; + template class SymmetricTensor<2, deal_II_dimension, deal_II_dimension+1>; +#endif +\} + +} + +for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS) { template void FEValuesViews::Scalar::get_function_values @@ -116,7 +164,7 @@ for (VEC : SERIAL_VECTORS) // Instantiations of functions in FEValuesBase -for (VEC : SERIAL_VECTORS) +for (VEC : SERIAL_VECTORS; deal_II_dimension : DIMENSIONS) { template void FEValuesBase::get_function_values @@ -341,3 +389,4 @@ for (VEC : SERIAL_VECTORS) #endif } + diff --git a/deal.II/source/fe/mapping.cc b/deal.II/source/fe/mapping.cc index e0d8290ede..97b1fba103 100644 --- a/deal.II/source/fe/mapping.cc +++ b/deal.II/source/fe/mapping.cc @@ -48,8 +48,8 @@ Mapping::transform_covariant ( const typename Mapping::InternalDataBase &mapping_data) const { Assert (offset == 0, ExcInternalError()); - - transform(input, output, mapping_data, mapping_covariant); + + transform(input, output, mapping_data, mapping_covariant); } @@ -116,10 +116,8 @@ Mapping::InternalDataBase::memory_consumption () const -template class Mapping; +// explicit instantiations +#include "mapping.inst" -#if deal_II_dimension != 3 -template class Mapping; -#endif DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping.inst.in b/deal.II/source/fe/mapping.inst.in new file mode 100644 index 0000000000..4417789b3f --- /dev/null +++ b/deal.II/source/fe/mapping.inst.in @@ -0,0 +1,23 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class Mapping; + +#if deal_II_dimension != 3 + template class Mapping; +#endif + } + diff --git a/deal.II/source/fe/mapping_c1.cc b/deal.II/source/fe/mapping_c1.cc index 215d63c221..d946b1394a 100644 --- a/deal.II/source/fe/mapping_c1.cc +++ b/deal.II/source/fe/mapping_c1.cc @@ -31,7 +31,6 @@ MappingC1::MappingC1 () } -#if deal_II_dimension == 1 template <> void @@ -42,11 +41,8 @@ MappingC1<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, Assert (dim > 1, ExcImpossibleInDim(dim)); } -#endif -#if deal_II_dimension == 2 - template <> void MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &cell, @@ -54,7 +50,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce { const unsigned int dim = 2; std::vector > line_points (2); - + // loop over each of the lines, // and if it is at the // boundary, then first get the @@ -152,7 +148,7 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce (face_vertex_normals[1][1] * std::cos(alpha) -face_vertex_normals[1][0] * std::sin(alpha))) -2*c; - + // next evaluate the so // determined cubic @@ -186,8 +182,6 @@ MappingC1<2>::add_line_support_points (const Triangulation<2>::cell_iterator &ce }; } -#endif - template @@ -200,9 +194,6 @@ MappingC1::add_line_support_points (const typename Triangulation void MappingC1<1>::add_quad_support_points (const Triangulation<1>::cell_iterator &, @@ -212,12 +203,8 @@ MappingC1<1>::add_quad_support_points (const Triangulation<1>::cell_iterator &, Assert (dim > 2, ExcImpossibleInDim(dim)); } -#endif - -#if deal_II_dimension == 2 - template <> void MappingC1<2>::add_quad_support_points (const Triangulation<2>::cell_iterator &, @@ -227,9 +214,6 @@ MappingC1<2>::add_quad_support_points (const Triangulation<2>::cell_iterator &, Assert (dim > 2, ExcImpossibleInDim(dim)); } -#endif - - template @@ -253,6 +237,7 @@ MappingC1::clone () const // explicit instantiations -template class MappingC1; +#include "mapping_c1.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_c1.inst.in b/deal.II/source/fe/mapping_c1.inst.in new file mode 100644 index 0000000000..27a107a0d9 --- /dev/null +++ b/deal.II/source/fe/mapping_c1.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class MappingC1; + } + diff --git a/deal.II/source/fe/mapping_cartesian.cc b/deal.II/source/fe/mapping_cartesian.cc index 39196105ab..7a06067feb 100644 --- a/deal.II/source/fe/mapping_cartesian.cc +++ b/deal.II/source/fe/mapping_cartesian.cc @@ -70,7 +70,7 @@ MappingCartesian::update_each (const UpdateFlags in) const UpdateFlags out = in; if (out & update_boundary_forms) out |= update_normal_vectors; - + return out; } @@ -86,7 +86,7 @@ MappingCartesian::get_data (const UpdateFlags update_flags, data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); data->update_flags = data->update_once | data->update_each; - + return data; } @@ -103,7 +103,7 @@ MappingCartesian::get_face_data (const UpdateFlags update_flags, data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); data->update_flags = data->update_once | data->update_each; - + return data; } @@ -120,7 +120,7 @@ MappingCartesian::get_subface_data (const UpdateFlags update_flag data->update_once = update_once(update_flags); data->update_each = update_each(update_flags); data->update_flags = data->update_once | data->update_each; - + return data; } @@ -157,7 +157,6 @@ MappingCartesian::compute_fill (const typename Triangulationface(face_no)->has_children() || (sub_no+1 < GeometryInfo::max_children_per_face+1), @@ -167,7 +166,6 @@ MappingCartesian::compute_fill (const typename Triangulationface(face_no)->has_children() || (sub_no < cell->face(face_no)->n_children()), ExcIndexRange (sub_no, 0, cell->face(face_no)->n_children())); -#endif } else // invalid face number, so @@ -180,7 +178,7 @@ MappingCartesian::compute_fill (const typename Triangulation start = cell->vertex(0); - + // Compute start point and sizes // along axes. Strange vertex // numbering makes this complicated @@ -205,7 +203,7 @@ MappingCartesian::compute_fill (const typename Triangulation::compute_fill (const typename Triangulation::faces_per_cell, ExcInternalError()); - + switch (dim) { case 2: @@ -302,7 +300,7 @@ template void MappingCartesian:: fill_fe_values (const typename Triangulation::cell_iterator& cell, - const Quadrature& q, + const Quadrature& q, typename Mapping::InternalDataBase& mapping_data, std::vector >& quadrature_points, std::vector& JxW_values, @@ -320,7 +318,7 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, InternalData &data = static_cast (mapping_data); std::vector > dummy; - + compute_fill (cell, invalid_face_number, invalid_face_number, cell_similarity, data, quadrature_points, @@ -359,8 +357,8 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, if (cell_similarity != CellSimilarity::translation) for (unsigned int i=0; i(); - // "compute" inverse Jacobian at the - // quadrature points, which are all + // "compute" inverse Jacobian at the + // quadrature points, which are all // the same if (data.current_update_flags() & update_inverse_jacobians) if (cell_similarity != CellSimilarity::translation) @@ -374,45 +372,77 @@ fill_fe_values (const typename Triangulation::cell_iterator& cell, -#if (deal_II_dimension == 1) +template <> +void +MappingCartesian<1,1>::fill_fe_face_values ( + const typename Triangulation<1,1>::cell_iterator &, + const unsigned int, + const Quadrature<0>&, + typename Mapping<1,1>::InternalDataBase&, + std::vector >&, + std::vector&, + std::vector >&, + std::vector >&) const +{ + Assert(false, ExcNotImplemented()); +} + -template +template <> void -MappingCartesian::fill_fe_face_values ( - const typename Triangulation::cell_iterator &, +MappingCartesian<1,2>::fill_fe_face_values ( + const typename Triangulation<1,2>::cell_iterator &, const unsigned int, - const Quadrature&, - typename Mapping::InternalDataBase&, - std::vector >&, + const Quadrature<0>&, + typename Mapping<1,2>::InternalDataBase&, + std::vector >&, std::vector&, - std::vector >&, - std::vector >&) const + std::vector >&, + std::vector >&) const { Assert(false, ExcNotImplemented()); } -template + +template <> +void +MappingCartesian<1,1>::fill_fe_subface_values ( + const typename Triangulation<1,1>::cell_iterator &, + const unsigned int, + const unsigned int, + const Quadrature<0>&, + typename Mapping<1,1>::InternalDataBase&, + std::vector >&, + std::vector&, + std::vector >&, + std::vector >&) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template <> void -MappingCartesian::fill_fe_subface_values ( - const typename Triangulation::cell_iterator &, +MappingCartesian<1,2>::fill_fe_subface_values ( + const typename Triangulation<1,2>::cell_iterator &, const unsigned int, const unsigned int, - const Quadrature&, - typename Mapping::InternalDataBase&, - std::vector >&, + const Quadrature<0>&, + typename Mapping<1,2>::InternalDataBase&, + std::vector >&, std::vector&, - std::vector >&, - std::vector >&) const + std::vector >&, + std::vector >&) const { Assert(false, ExcNotImplemented()); } -#else // Implementation for dim != 1 - + template void MappingCartesian::fill_fe_face_values ( @@ -433,7 +463,7 @@ MappingCartesian::fill_fe_face_values ( ExcInternalError()); InternalData &data = static_cast (mapping_data); - compute_fill (cell, face_no, invalid_face_number, + compute_fill (cell, face_no, invalid_face_number, CellSimilarity::none, data, quadrature_points, @@ -447,7 +477,7 @@ MappingCartesian::fill_fe_face_values ( for (unsigned int d=0; d::unit_normal_direction[face_no]) J *= data.length[d]; - + if (data.current_update_flags() & update_JxW_values) for (unsigned int i=0; i::fill_fe_subface_values ( for (unsigned int d=0; d::unit_normal_direction[face_no]) J *= data.length[d]; - + if (data.current_update_flags() & update_JxW_values) { // Here, @@ -534,10 +564,8 @@ MappingCartesian::fill_fe_subface_values ( } } - -#endif - + template void MappingCartesian::transform ( @@ -547,28 +575,28 @@ MappingCartesian::transform ( const MappingType mapping_type) const { AssertDimension (input.size(), output.size()); - Assert (dynamic_cast(&mapping_data) != 0, - ExcInternalError()); + Assert (dynamic_cast(&mapping_data) != 0, + ExcInternalError()); const InternalData &data = static_cast (mapping_data); - + switch (mapping_type) { case mapping_covariant: { Assert (data.update_flags & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( typename FEValuesBase::ExcAccessToUninitializedField()); Assert (data.update_flags & update_volume_elements, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( } } + + template void MappingCartesian::transform ( @@ -600,29 +630,29 @@ MappingCartesian::transform ( const MappingType mapping_type) const { AssertDimension (input.size(), output.size()); - Assert (dynamic_cast(&mapping_data) != 0, + Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast(mapping_data); - + switch (mapping_type) { case mapping_covariant: { Assert (data.update_flags & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( { Assert (data.update_flags & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( { Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( typename FEValuesBase::ExcAccessToUninitializedField()); Assert (data.update_flags & update_volume_elements, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( / data.length[d1] / data.volume_element; return; } - + default: Assert(false, ExcNotImplemented()); } @@ -755,7 +785,7 @@ MappingCartesian::clone () const //--------------------------------------------------------------------------- // explicit instantiations +#include "mapping_cartesian.inst" -template class MappingCartesian; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_cartesian.inst.in b/deal.II/source/fe/mapping_cartesian.inst.in new file mode 100644 index 0000000000..ab9c3cf951 --- /dev/null +++ b/deal.II/source/fe/mapping_cartesian.inst.in @@ -0,0 +1,19 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class MappingCartesian; + } + diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index 33253296db..e4e54d0288 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -56,8 +56,6 @@ MappingQ::InternalData::memory_consumption () const -#if deal_II_dimension == 1 - // in 1d, it is irrelevant which polynomial degree to use, since all // cells are scaled linearly. Unless codimension is equal to two template<> @@ -92,8 +90,6 @@ template<> MappingQ<1>::~MappingQ () {} -#endif - namespace @@ -182,8 +178,6 @@ MappingQ::~MappingQ () -#if deal_II_dimension == 1 - template<> void MappingQ<1>::compute_shapes_virtual (const std::vector > &unit_points, @@ -192,8 +186,6 @@ MappingQ<1>::compute_shapes_virtual (const std::vector > &unit_points, MappingQ1<1>::compute_shapes_virtual(unit_points, data); } -#endif - template @@ -489,7 +481,6 @@ MappingQ::fill_fe_subface_values (const typename Triangulation void @@ -498,7 +489,7 @@ MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const Assert(false, ExcInternalError()); } -#else + template void @@ -704,10 +695,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const ExcInternalError()); } -#endif - -#if deal_II_dimension==3 template <> void @@ -752,7 +740,6 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const ExcInternalError()); } -#endif template @@ -764,9 +751,6 @@ MappingQ::set_laplace_on_hex_vector(Table<2,double> &) const - -#if deal_II_dimension==1 - template <> void MappingQ<1>::compute_laplace_vector(Table<2,double> &) const @@ -774,7 +758,6 @@ MappingQ<1>::compute_laplace_vector(Table<2,double> &) const Assert(false, ExcInternalError()); } -#else template @@ -836,8 +819,6 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const lvs(i,k) = -S_1_T(i,k); } -#endif - template @@ -959,10 +940,6 @@ MappingQ::compute_support_points_laplace(const typename Triangulat - - -#if deal_II_dimension==1 - template <> void MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, @@ -975,6 +952,7 @@ MappingQ<1>::add_line_support_points (const Triangulation<1>::cell_iterator &, } + template <> void MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator &cell, @@ -1006,7 +984,6 @@ MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator } } -#endif template @@ -1080,10 +1057,6 @@ MappingQ::add_line_support_points (const typename Triangulation void MappingQ<3>:: @@ -1250,9 +1223,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, } } -#endif -#if deal_II_dimension == 2 template<> void @@ -1268,7 +1239,6 @@ add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, a.push_back(quad_points[i]); } -#endif template @@ -1281,6 +1251,7 @@ add_quad_support_points(const typename Triangulation::cell_iterato } + template void MappingQ::transform ( @@ -1319,6 +1290,7 @@ MappingQ::transform ( } + template void MappingQ::transform ( @@ -1443,11 +1415,9 @@ MappingQ::clone () const } -// explicit instantiation -template class MappingQ; -#if deal_II_dimension < 3 -template class MappingQ; -#endif +// explicit instantiations +#include "mapping_q.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_q.inst.in b/deal.II/source/fe/mapping_q.inst.in new file mode 100644 index 0000000000..d21499b464 --- /dev/null +++ b/deal.II/source/fe/mapping_q.inst.in @@ -0,0 +1,23 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class MappingQ; + +#if deal_II_dimension < 3 + template class MappingQ; +#endif + } + diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index ee72c8f37b..2b66511b70 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -863,7 +863,6 @@ MappingQ1::fill_fe_values ( } -#if deal_II_dimension == 2 template<> void @@ -882,64 +881,117 @@ MappingQ1<2,3>::compute_fill_face (const Triangulation<2,3>::cell_iterator &, Assert(false, ExcNotImplemented()); } -#endif -#if (deal_II_dimension == 1) -template +template <> void -MappingQ1::compute_fill_face ( - const typename Triangulation::cell_iterator &, +MappingQ1<1,1>::compute_fill_face ( + const typename Triangulation<1,1>::cell_iterator &, + const unsigned int, + const unsigned int, + const unsigned int, + const DataSetDescriptor, + const std::vector &, + InternalData &, + std::vector > &, + std::vector &, + std::vector > &, + std::vector > &) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template <> +void +MappingQ1<1,2>::compute_fill_face ( + const typename Triangulation<1,2>::cell_iterator &, const unsigned int, const unsigned int, const unsigned int, const DataSetDescriptor, const std::vector &, InternalData &, - std::vector > &, + std::vector > &, std::vector &, - std::vector > &, - std::vector > &) const + std::vector > &, + std::vector > &) const { Assert(false, ExcNotImplemented()); } -template + +template <> void -MappingQ1::fill_fe_face_values ( - const typename Triangulation::cell_iterator &, +MappingQ1<1,1>::fill_fe_face_values ( + const typename Triangulation<1,1>::cell_iterator &, const unsigned, - const Quadrature&, - typename Mapping::InternalDataBase&, - std::vector >&, + const Quadrature<0>&, + typename Mapping<1,1>::InternalDataBase&, + std::vector >&, std::vector&, - std::vector >&, - std::vector >&) const + std::vector >&, + std::vector >&) const { Assert(false, ExcNotImplemented()); } -template + +template <> void -MappingQ1::fill_fe_subface_values ( - const typename Triangulation::cell_iterator &, +MappingQ1<1,2>::fill_fe_face_values ( + const typename Triangulation<1,2>::cell_iterator &, + const unsigned, + const Quadrature<0>&, + typename Mapping<1,2>::InternalDataBase&, + std::vector >&, + std::vector&, + std::vector >&, + std::vector >&) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template <> +void +MappingQ1<1,1>::fill_fe_subface_values ( + const typename Triangulation<1,1>::cell_iterator &, const unsigned, const unsigned, - const Quadrature&, - typename Mapping::InternalDataBase&, - std::vector >&, + const Quadrature<0>&, + typename Mapping<1,1>::InternalDataBase&, + std::vector >&, std::vector&, - std::vector >&, - std::vector >&) const + std::vector >&, + std::vector >&) const +{ + Assert(false, ExcNotImplemented()); +} + + + +template <> +void +MappingQ1<1,2>::fill_fe_subface_values ( + const typename Triangulation<1,2>::cell_iterator &, + const unsigned, + const unsigned, + const Quadrature<0>&, + typename Mapping<1,2>::InternalDataBase&, + std::vector >&, + std::vector&, + std::vector >&, + std::vector >&) const { Assert(false, ExcNotImplemented()); } -#else -// Implementation for 2D and 3D template void @@ -1102,8 +1154,6 @@ MappingQ1::fill_fe_subface_values (const typename Triangulation void @@ -1369,7 +1419,6 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it } -#if deal_II_dimension == 1 template<> void @@ -1383,9 +1432,7 @@ transform_real_to_unit_cell_internal Assert(false, ExcNotImplemented()); } -#endif -#if deal_II_dimension == 2 template<> void @@ -1399,7 +1446,6 @@ transform_real_to_unit_cell_internal Assert(false, ExcNotImplemented()); } -#endif template @@ -1486,13 +1532,8 @@ MappingQ1 StaticMappingQ1::mapping; -//--------------------------------------------------------------------------- +//--------------------------- Explicit instantiations ----------------------- +#include "mapping_q1.inst" -template class MappingQ1; -template struct StaticMappingQ1; -#if deal_II_dimension != 3 -template class MappingQ1; -template struct StaticMappingQ1; -#endif DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_q1.inst.in b/deal.II/source/fe/mapping_q1.inst.in new file mode 100644 index 0000000000..04e52551ae --- /dev/null +++ b/deal.II/source/fe/mapping_q1.inst.in @@ -0,0 +1,25 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class MappingQ1; + template struct StaticMappingQ1; + +#if deal_II_dimension != 3 + template class MappingQ1; + template struct StaticMappingQ1; +#endif + } + diff --git a/deal.II/source/fe/mapping_q1_eulerian.cc b/deal.II/source/fe/mapping_q1_eulerian.cc index 3dd1292162..4c67d46f39 100644 --- a/deal.II/source/fe/mapping_q1_eulerian.cc +++ b/deal.II/source/fe/mapping_q1_eulerian.cc @@ -139,19 +139,8 @@ MappingQ1Eulerian::fill_fe_values ( -// explicit instantiation -template class MappingQ1Eulerian >; -template class MappingQ1Eulerian >; -#ifdef DEAL_II_USE_PETSC -template class MappingQ1Eulerian; -#endif - -// Explicit instantiation for codimension one problems. -#if deal_II_dimension != 3 -template class MappingQ1Eulerian, deal_II_dimension+1 >; -# ifdef DEAL_II_USE_PETSC -template class MappingQ1Eulerian; -# endif -#endif +// explicit instantiations +#include "mapping_q1_eulerian.inst" + DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_q1_eulerian.inst.in b/deal.II/source/fe/mapping_q1_eulerian.inst.in new file mode 100644 index 0000000000..3c822d119b --- /dev/null +++ b/deal.II/source/fe/mapping_q1_eulerian.inst.in @@ -0,0 +1,30 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class MappingQ1Eulerian >; + template class MappingQ1Eulerian >; +#ifdef DEAL_II_USE_PETSC + template class MappingQ1Eulerian; +#endif + +#if deal_II_dimension != 3 + template class MappingQ1Eulerian, deal_II_dimension+1 >; +# ifdef DEAL_II_USE_PETSC + template class MappingQ1Eulerian; +# endif +#endif + } + diff --git a/deal.II/source/fe/mapping_q_eulerian.cc b/deal.II/source/fe/mapping_q_eulerian.cc index f0ded0c90b..6a1e9010cd 100644 --- a/deal.II/source/fe/mapping_q_eulerian.cc +++ b/deal.II/source/fe/mapping_q_eulerian.cc @@ -208,22 +208,8 @@ MappingQEulerian::fill_fe_values ( -// explicit instantiation -template class MappingQEulerian >; -#ifdef DEAL_II_USE_PETSC -template class MappingQEulerian; -#endif - -#if deal_II_dimension != 3 -template class MappingQEulerian, - deal_II_dimension+1>; - -# ifdef DEAL_II_USE_PETSC -template class MappingQEulerian; -# endif - -#endif +// explicit instantiations +#include "mapping_q_eulerian.inst" DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/fe/mapping_q_eulerian.inst.in b/deal.II/source/fe/mapping_q_eulerian.inst.in new file mode 100644 index 0000000000..b8e323ca3d --- /dev/null +++ b/deal.II/source/fe/mapping_q_eulerian.inst.in @@ -0,0 +1,33 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 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. +// +//--------------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS) + { + template class MappingQEulerian >; +#ifdef DEAL_II_USE_PETSC + template class MappingQEulerian; +#endif + +#if deal_II_dimension != 3 + template class MappingQEulerian, + deal_II_dimension+1>; + +# ifdef DEAL_II_USE_PETSC + template class MappingQEulerian; +# endif + +#endif + } + -- 2.39.5