From f6a17c2ec3e9086fee0de1e054393bccb4ffec99 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 3 Nov 2010 03:34:27 +0000 Subject: [PATCH] Re-apply previous. It now compiles, links and should even work with some luck. git-svn-id: https://svn.dealii.org/trunk@22592 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 + deal.II/include/deal.II/fe/mapping.h | 56 +++--- deal.II/include/deal.II/fe/mapping_q.h | 8 +- deal.II/include/deal.II/fe/mapping_q1.h | 12 +- deal.II/include/deal.II/hp/fe_values.h | 36 ++-- .../deal.II/numerics/vectors.templates.h | 23 ++- deal.II/source/fe/fe_values.cc | 20 +- deal.II/source/fe/mapping_q.cc | 140 ++++++------- deal.II/source/fe/mapping_q1.cc | 189 +++++++++--------- deal.II/source/hp/fe_values.cc | 82 ++++---- deal.II/source/numerics/vectors.cc | 130 ++++++++++-- 11 files changed, 404 insertions(+), 301 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index ea5e92646e..c6e2146b58 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -332,6 +332,15 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.

deal.II

    +

  1. New: The class hp::FEFaceValues and hp::FESubfaceValues were not + previously available for meshes that were embedded in a higher + dimensional space (i.e. if the codimension was greater than 1). This is + now fixed. As a consequence, the VectorTools::interpolate_boundary_values + function is now also available for such meshes. +
    + (WB, 2010/11/02) +

    2. +

  2. Fixed: The FEValuesExtractors::Vector class did not work when the dimension of the domain was not equal to the dimension of the space in which it is embedded. This is now fixed. diff --git a/deal.II/include/deal.II/fe/mapping.h b/deal.II/include/deal.II/fe/mapping.h index 4cc31d3ecc..d3dc65cb74 100644 --- a/deal.II/include/deal.II/fe/mapping.h +++ b/deal.II/include/deal.II/fe/mapping.h @@ -126,7 +126,7 @@ enum MappingType * * A hint to implementators: no function except the two functions * @p update_once and @p update_each may add any flags. - * + * * For more information about the spacedim template parameter * check the documentation of FiniteElement or the one of * Triangulation. @@ -138,12 +138,12 @@ template class Mapping : public Subscriptor { public: - + /** * Virtual destructor. */ virtual ~Mapping (); - + /** * Transforms the point @p p on * the unit cell to the point @@ -154,7 +154,7 @@ class Mapping : public Subscriptor transform_unit_to_real_cell ( const typename Triangulation::cell_iterator &cell, const Point &p) const = 0; - + /** * Transforms the point @p p on * the real cell to the point @@ -165,7 +165,7 @@ class Mapping : public Subscriptor transform_real_to_unit_cell ( const typename Triangulation::cell_iterator &cell, const Point &p) const = 0; - + /** * Base class for internal data * of finite element and mapping @@ -219,7 +219,7 @@ class Mapping : public Subscriptor * @p first_cell to @p true. */ InternalDataBase (); - + /** * Virtual destructor for * derived classes @@ -231,7 +231,7 @@ class Mapping : public Subscriptor * by reinit. */ UpdateFlags update_flags; - + /** * Values computed by * constructor. @@ -286,7 +286,7 @@ class Mapping : public Subscriptor * the first cell. */ virtual void clear_first_cell (); - + /** * Return an estimate (in * bytes) or the memory @@ -302,14 +302,14 @@ class Mapping : public Subscriptor * if #update_volume_elements. */ std::vector volume_elements; - + /** * The positions of the * mapped (generalized) * support points. */ std::vector > support_point_values; - + /* * The Jacobian of the * transformation in the @@ -317,7 +317,7 @@ class Mapping : public Subscriptor * points. */ std::vector > support_point_gradients; - + /* * The inverse of the * Jacobian of the @@ -326,8 +326,8 @@ class Mapping : public Subscriptor * points. */ std::vector > support_point_inverse_gradients; - - + + private: /** * The value returned by @@ -335,7 +335,7 @@ class Mapping : public Subscriptor */ bool first_cell; }; - + /** * Transform a field of * vectors accorsing to @@ -360,7 +360,7 @@ class Mapping : public Subscriptor VectorSlice > > output, const InternalDataBase &internal, const MappingType type) const = 0; - + /** * @deprecated Use transform() instead. */ @@ -391,7 +391,7 @@ class Mapping : public Subscriptor /** * @deprecated Use transform() instead. */ - + void transform_contravariant (const VectorSlice > > intput, const unsigned int offset, @@ -405,7 +405,7 @@ class Mapping : public Subscriptor const Point& support_point_value( const unsigned int index, const typename Mapping::InternalDataBase &internal) const; - + /** * The Jacobian * matrix of the transformation @@ -415,7 +415,7 @@ class Mapping : public Subscriptor const Tensor<2,spacedim>& support_point_gradient( const unsigned int index, const typename Mapping::InternalDataBase &internal) const; - + /** * The inverse Jacobian * matrix of the transformation @@ -425,7 +425,7 @@ class Mapping : public Subscriptor const Tensor<2,spacedim>& support_point_inverse_gradient( const unsigned int index, const typename Mapping::InternalDataBase &internal) const; - + /** * Return a pointer to a copy of the * present object. The caller of this @@ -441,7 +441,7 @@ class Mapping : public Subscriptor */ virtual Mapping * clone () const = 0; - + /** * Returns whether the mapping preserves * vertex locations, i.e. whether the @@ -466,7 +466,7 @@ class Mapping : public Subscriptor DeclException0 (ExcInvalidData); private: - + /** * Indicate fields to be updated * in the constructor of @@ -479,7 +479,7 @@ class Mapping : public Subscriptor * See @ref UpdateFlagsEssay. */ virtual UpdateFlags update_once (const UpdateFlags) const = 0; - + /** * The same as update_once(), * but for the flags to be updated for @@ -507,7 +507,7 @@ class Mapping : public Subscriptor virtual InternalDataBase* get_face_data (const UpdateFlags flags, const Quadrature& quadrature) const = 0; - + /** * Prepare internal data * structure for transformation @@ -601,9 +601,9 @@ class Mapping : public Subscriptor const unsigned int face_no, const Quadrature &quadrature, InternalDataBase &internal, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &boundary_form, + std::vector > &boundary_form, std::vector > &normal_vectors) const = 0; /** @@ -615,9 +615,9 @@ class Mapping : public Subscriptor const unsigned int sub_no, const Quadrature &quadrature, InternalDataBase &internal, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &boundary_form, + std::vector > &boundary_form, std::vector > &normal_vectors) const = 0; /** @@ -682,7 +682,7 @@ Mapping::support_point_value( const typename Mapping::InternalDataBase& internal) const { AssertIndexRange(index, internal.support_point_values.size()); - return internal.support_point_values[index]; + return internal.support_point_values[index]; } diff --git a/deal.II/include/deal.II/fe/mapping_q.h b/deal.II/include/deal.II/fe/mapping_q.h index 4b063526b3..75538439c1 100644 --- a/deal.II/include/deal.II/fe/mapping_q.h +++ b/deal.II/include/deal.II/fe/mapping_q.h @@ -222,9 +222,9 @@ class MappingQ : public MappingQ1 const unsigned int face_no, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, + typename std::vector > &quadrature_points, std::vector &JxW_values, - typename std::vector > &exterior_form, + typename std::vector > &exterior_form, typename std::vector > &normal_vectors) const ; /** @@ -237,9 +237,9 @@ class MappingQ : public MappingQ1 const unsigned int sub_no, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, + typename std::vector > &quadrature_points, std::vector &JxW_values, - typename std::vector > &exterior_form, + typename std::vector > &exterior_form, typename std::vector > &normal_vectors) const ; /** diff --git a/deal.II/include/deal.II/fe/mapping_q1.h b/deal.II/include/deal.II/fe/mapping_q1.h index 616d535ade..33ff0e7e7e 100644 --- a/deal.II/include/deal.II/fe/mapping_q1.h +++ b/deal.II/include/deal.II/fe/mapping_q1.h @@ -312,9 +312,9 @@ class MappingQ1 : public Mapping const unsigned int face_no, const Quadrature &quadrature, typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, + typename std::vector > &quadrature_points, std::vector &JxW_values, - typename std::vector > &boundary_form, + typename std::vector > &boundary_form, typename std::vector > &normal_vectors) const ; /** @@ -327,9 +327,9 @@ class MappingQ1 : public Mapping const unsigned int sub_no, const Quadrature& quadrature, typename Mapping::InternalDataBase &mapping_data, - typename std::vector > &quadrature_points, + typename std::vector > &quadrature_points, std::vector &JxW_values, - typename std::vector > &boundary_form, + typename std::vector > &boundary_form, typename std::vector > &normal_vectors) const ; /** @@ -399,9 +399,9 @@ class MappingQ1 : public Mapping const DataSetDescriptor data_set, const std::vector &weights, InternalData &mapping_data, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &boundary_form, + std::vector > &boundary_form, std::vector > &normal_vectors) const; /** diff --git a/deal.II/include/deal.II/hp/fe_values.h b/deal.II/include/deal.II/hp/fe_values.h index f39cbc0874..308f9946b0 100644 --- a/deal.II/include/deal.II/hp/fe_values.h +++ b/deal.II/include/deal.II/hp/fe_values.h @@ -546,8 +546,8 @@ namespace hp * @ingroup hp hpcollection * @author Wolfgang Bangerth, 2003 */ - template - class FEFaceValues : public internal::hp::FEValuesBase > + template + class FEFaceValues : public internal::hp::FEValuesBase > { public: /** @@ -570,8 +570,8 @@ namespace hp * DoFHandler::get_fe() * function. */ - FEFaceValues (const hp::MappingCollection &mapping_collection, - const hp::FECollection &fe_collection, + FEFaceValues (const hp::MappingCollection &mapping_collection, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags); @@ -598,7 +598,7 @@ namespace hp * DoFHandler::get_fe() * function. */ - FEFaceValues (const hp::FECollection &fe_collection, + FEFaceValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags); @@ -702,7 +702,7 @@ namespace hp * this argument is specified. */ void - reinit (const typename hp::DoFHandler::cell_iterator &cell, + reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index = numbers::invalid_unsigned_int, const unsigned int mapping_index = numbers::invalid_unsigned_int, @@ -738,7 +738,7 @@ namespace hp * these last three arguments. */ void - reinit (const typename dealii::DoFHandler::cell_iterator &cell, + reinit (const typename dealii::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index = numbers::invalid_unsigned_int, const unsigned int mapping_index = numbers::invalid_unsigned_int, @@ -774,7 +774,7 @@ namespace hp * these last three arguments. */ void - reinit (const typename MGDoFHandler::cell_iterator &cell, + reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index = numbers::invalid_unsigned_int, const unsigned int mapping_index = numbers::invalid_unsigned_int, @@ -811,7 +811,7 @@ namespace hp * three arguments. */ void - reinit (const typename Triangulation::cell_iterator &cell, + reinit (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index = numbers::invalid_unsigned_int, const unsigned int mapping_index = numbers::invalid_unsigned_int, @@ -827,8 +827,8 @@ namespace hp * @ingroup hp hpcollection * @author Wolfgang Bangerth, 2003 */ - template - class FESubfaceValues : public internal::hp::FEValuesBase > + template + class FESubfaceValues : public internal::hp::FEValuesBase > { public: /** @@ -851,8 +851,8 @@ namespace hp * DoFHandler::get_fe() * function. */ - FESubfaceValues (const hp::MappingCollection &mapping_collection, - const hp::FECollection &fe_collection, + FESubfaceValues (const hp::MappingCollection &mapping_collection, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags); @@ -879,7 +879,7 @@ namespace hp * DoFHandler::get_fe() * function. */ - FESubfaceValues (const hp::FECollection &fe_collection, + FESubfaceValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags); @@ -962,7 +962,7 @@ namespace hp * this argument is specified. */ void - reinit (const typename hp::DoFHandler::cell_iterator &cell, + reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index = numbers::invalid_unsigned_int, @@ -999,7 +999,7 @@ namespace hp * these last three arguments. */ void - reinit (const typename dealii::DoFHandler::cell_iterator &cell, + reinit (const typename dealii::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index = numbers::invalid_unsigned_int, @@ -1036,7 +1036,7 @@ namespace hp * these last three arguments. */ void - reinit (const typename MGDoFHandler::cell_iterator &cell, + reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index = numbers::invalid_unsigned_int, @@ -1074,7 +1074,7 @@ namespace hp * three arguments. */ void - reinit (const typename Triangulation::cell_iterator &cell, + reinit (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index = numbers::invalid_unsigned_int, diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h index e0494237be..b27f75adf4 100644 --- a/deal.II/include/deal.II/numerics/vectors.templates.h +++ b/deal.II/include/deal.II/numerics/vectors.templates.h @@ -1353,6 +1353,7 @@ interpolate_boundary_values (const Mapping const std::vector &component_mask_) { const unsigned int dim=DH::dimension; + const unsigned int spacedim=DH::space_dimension; Assert ((component_mask_.size() == 0) || (component_mask_.size() == dof.get_fe().n_components()), @@ -1373,7 +1374,7 @@ interpolate_boundary_values (const Mapping const unsigned int n_components = DoFTools::n_components(dof); const bool fe_is_system = (n_components != 1); - for (typename FunctionMap::type::const_iterator i=function_map.begin(); + for (typename FunctionMap::type::const_iterator i=function_map.begin(); i!=function_map.end(); ++i) Assert (n_components == i->second->n_components, ExcDimensionMismatch(n_components, i->second->n_components)); @@ -1391,7 +1392,7 @@ interpolate_boundary_values (const Mapping std::vector face_dofs; face_dofs.reserve (DoFTools::max_dofs_per_face(dof)); - std::vector > dof_locations; + std::vector > dof_locations; dof_locations.reserve (DoFTools::max_dofs_per_face(dof)); // array to store the values of @@ -1411,11 +1412,11 @@ interpolate_boundary_values (const Mapping // the interpolation points of all // finite elements that may ever be // in use - hp::FECollection finite_elements (dof.get_fe()); + hp::FECollection finite_elements (dof.get_fe()); hp::QCollection q_collection; for (unsigned int f=0; f &fe = finite_elements[f]; + const FiniteElement &fe = finite_elements[f]; // generate a quadrature rule // on the face from the unit @@ -1483,9 +1484,9 @@ interpolate_boundary_values (const Mapping // hp::FEFaceValues object that we // can use to evaluate the boundary // values at - hp::MappingCollection mapping_collection (mapping); - hp::FEFaceValues x_fe_values (mapping_collection, finite_elements, q_collection, - update_quadrature_points); + hp::MappingCollection mapping_collection (mapping); + hp::FEFaceValues x_fe_values (mapping_collection, finite_elements, q_collection, + update_quadrature_points); typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); @@ -1494,7 +1495,7 @@ interpolate_boundary_values (const Mapping for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; ++face_no) { - const FiniteElement &fe = cell->get_fe(); + const FiniteElement &fe = cell->get_fe(); // we can presently deal only with // primitive elements for boundary @@ -1526,14 +1527,14 @@ interpolate_boundary_values (const Mapping { // face is of the right component x_fe_values.reinit(cell, face_no); - const FEFaceValues &fe_values = x_fe_values.get_present_fe_values(); + const FEFaceValues &fe_values = x_fe_values.get_present_fe_values(); // get indices, physical location and // boundary values of dofs on this // face face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); - const std::vector > &dof_locations + const std::vector > &dof_locations = fe_values.get_quadrature_points (); if (fe_is_system) @@ -1997,7 +1998,7 @@ VectorTools::project_boundary_values (const Mapping &mappin && ! excluded_dofs[dof_to_boundary_mapping[i]]) { Assert(numbers::is_finite(boundary_projection(dof_to_boundary_mapping[i])), ExcNumberNotFinite()); - + // this dof is on one of the // interesting boundary parts // diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc index e89b8a4911..8481733c2f 100644 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@ -3607,7 +3607,7 @@ FEFaceValues::FEFaceValues (const FiniteElement &fe, FEFaceValuesBase (quadrature.size(), fe.dofs_per_cell, update_flags, - StaticMappingQ1::mapping, + StaticMappingQ1::mapping, fe, quadrature) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -3850,7 +3850,7 @@ FESubfaceValues::FESubfaceValues (const FiniteElement (quadrature.size(), fe.dofs_per_cell, update_flags, - StaticMappingQ1::mapping, + StaticMappingQ1::mapping, fe, quadrature) { Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping")); @@ -4158,14 +4158,14 @@ 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 +template class FEValuesBase:: + CellIterator::cell_iterator>; + +#if deal_II_dimension == 2 +template class FEFaceValuesBase; +template class FEFaceValues; +template class FESubfaceValues; +#endif #endif diff --git a/deal.II/source/fe/mapping_q.cc b/deal.II/source/fe/mapping_q.cc index 82bc61a809..33253296db 100644 --- a/deal.II/source/fe/mapping_q.cc +++ b/deal.II/source/fe/mapping_q.cc @@ -46,7 +46,7 @@ MappingQ::InternalData::InternalData (const unsigned int n_shape_f template unsigned int -MappingQ::InternalData::memory_consumption () const +MappingQ::InternalData::memory_consumption () const { return (MappingQ1::InternalData::memory_consumption () + MemoryConsumption::memory_consumption (unit_normals) + @@ -143,7 +143,7 @@ MappingQ::MappingQ (const unsigned int p, Assert (n_shape_functions==tensor_pols->n(), ExcInternalError()); Assert(n_inner+n_outer==n_shape_functions, ExcInternalError()); - + // build laplace_on_quad_vector if (degree>1) { @@ -216,7 +216,7 @@ MappingQ::compute_shapes_virtual (const std::vector > & ExcInternalError()); grads.resize(n_shape_functions); } - + // // dummy variable of size 0 std::vector > grad2; if (data.shape_second_derivatives.size()!=0) @@ -226,16 +226,16 @@ MappingQ::compute_shapes_virtual (const std::vector > & grad2.resize(n_shape_functions); } - + if (data.shape_values.size()!=0 || data.shape_derivatives.size()!=0) for (unsigned int point=0; pointcompute(unit_points[point], values, grads, grad2); - + if (data.shape_values.size()!=0) for (unsigned int i=0; i::fill_fe_values ( if (get_degree() > 1) cell_similarity = CellSimilarity::invalid_next_cell; } - + MappingQ1::fill_fe_values(cell, q, *p_data, quadrature_points, JxW_values, jacobians, jacobian_grads, inverse_jacobians, @@ -371,9 +371,9 @@ MappingQ::fill_fe_face_values ( const unsigned int face_no, const Quadrature &q, typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &exterior_forms, + std::vector > &exterior_forms, std::vector > &normal_vectors) const { // convert data object to internal @@ -383,7 +383,7 @@ MappingQ::fill_fe_face_values ( Assert (dynamic_cast (&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast (mapping_data); - + // check whether this cell needs // the full mapping or can be // treated by a reduced Q1 mapping, @@ -433,9 +433,9 @@ MappingQ::fill_fe_subface_values (const typename Triangulation &q, typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &exterior_forms, + std::vector > &exterior_forms, std::vector > &normal_vectors) const { // convert data object to internal @@ -518,7 +518,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const // for degree==1, we shouldn't have to // compute any support points, since // all of them are on the vertices - + case 2: { // (checked these values against the @@ -533,7 +533,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const Assert (sizeof(loqv2)/sizeof(loqv2[0]) == n_inner_2d * n_outer_2d, ExcInternalError()); - + break; } @@ -542,7 +542,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const // (same as above) static const double loqv3[4*12] ={80/1053., 1/81., 1/81., 11/1053., 25/117., 44/351., - 7/117., 16/351., 25/117., 44/351., 7/117., 16/351., + 7/117., 16/351., 25/117., 44/351., 7/117., 16/351., 1/81., 80/1053., 11/1053., 1/81., 7/117., 16/351., 25/117., 44/351., 44/351., 25/117., 16/351., 7/117., 1/81., 11/1053., 80/1053., 1/81., 44/351., 25/117., @@ -552,9 +552,9 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const Assert (sizeof(loqv3)/sizeof(loqv3[0]) == n_inner_2d * n_outer_2d, ExcInternalError()); - + loqv_ptr=&loqv3[0]; - + break; } @@ -570,7 +570,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.2231273865431891, 0.1346851306015187, 0.03812914216116723, 0.02913160002633253, 0.02200737428129391, 0.01600835564431222, - + 0.00664803151334206, 0.006648031513342719, 0.002873452861657458, 0.002873452861657626, 0.07903572682584378, 0.05969238281250031, @@ -588,7 +588,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.03812914216116729, 0.1346851306015185, 0.2231273865431898, 0.01600835564431217, 0.02200737428129394, 0.02913160002633262, - + 0.006648031513342073, 0.002873452861657473, 0.006648031513342726, 0.002873452861657636, 0.1527716818820238, 0.2348152760709273, @@ -597,7 +597,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.07903572682584376, 0.05969238281250026, 0.03619864817415824, 0.07903572682584187, 0.0596923828124998, 0.0361986481741581, - + 0.01106770833333302, 0.01106770833333336, 0.01106770833333337, 0.01106770833333374, 0.06770833333333424, 0.1035156250000011, @@ -606,7 +606,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.06770833333333422, 0.1035156250000009, 0.06770833333333436, 0.0677083333333337, 0.1035156249999988, 0.0677083333333339, - + 0.002873452861657185, 0.006648031513342362, 0.002873452861657334, 0.006648031513343038, 0.02496269311797779, 0.04081948955407401, @@ -615,7 +615,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.03619864817415819, 0.05969238281250028, 0.07903572682584407, 0.03619864817415804, 0.05969238281249986, 0.0790357268258422, - + -0.001075744628906913, 0.00191429223907134, 0.07405921850311592, -0.001075744628905865, 0.03812914216116729, 0.1346851306015185, @@ -633,7 +633,7 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.02496269311797776, 0.04081948955407392, 0.02496269311797785, 0.1527716818820233, 0.2348152760709258, 0.1527716818820236, - + 0.001914292239071237, -0.001075744628906803, -0.001075744628906778, 0.07405921850311617, 0.01600835564431228, 0.02200737428129401, @@ -641,24 +641,24 @@ MappingQ::set_laplace_on_quad_vector(Table<2,double> &loqvs) const 0.1346851306015182, 0.2231273865431886, 0.01600835564431228, 0.02200737428129397, 0.02913160002633523, 0.03812914216116726, - 0.1346851306015181, 0.2231273865431886, + 0.1346851306015181, 0.2231273865431886, }; - + Assert (sizeof(loqv4)/sizeof(loqv4[0]) == n_inner_2d * n_outer_2d, ExcInternalError()); - + loqv_ptr=&loqv4[0]; - + break; } - + // no other cases implemented, // so simply fall through default: break; } - + if (loqv_ptr!=0) { // precomputed. copy values to @@ -726,10 +726,10 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 7/192., 1/12., 1/12., 1/12., 1/12., 1/12., 1/12.}; - + lohv_ptr=&loqv2[0]; } - + if (lohv_ptr!=0) { // precomputed. copy values to @@ -742,7 +742,7 @@ MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const else // not precomputed, then do so now compute_laplace_vector(lohvs); - + // the sum of weights of the points // at the outer rim should be // one. check this @@ -794,11 +794,11 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const // points on the unit cell const QGauss quadrature(degree+1); const unsigned int n_q_points=quadrature.size(); - + InternalData quadrature_data(n_shape_functions); quadrature_data.shape_derivatives.resize(n_shape_functions * n_q_points); this->compute_shapes(quadrature.get_points(), quadrature_data); - + // Compute the stiffness matrix of // the inner dofs FullMatrix S(n_inner); @@ -808,7 +808,7 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const S(i,j) += contract(quadrature_data.derivative(point, n_outer+i), quadrature_data.derivative(point, n_outer+j)) * quadrature.weight(point); - + // Compute the components of T to be the // product of gradients of inner and // outer shape functions. @@ -819,15 +819,15 @@ MappingQ::compute_laplace_vector(Table<2,double> &lvs) const T(i,k) += contract(quadrature_data.derivative(point, n_outer+i), quadrature_data.derivative(point, k)) *quadrature.weight(point); - + FullMatrix S_1(n_inner); S_1.invert(S); - + FullMatrix S_1_T(n_inner, n_outer); - + // S:=S_1*T S_1.mmult(S_1_T,T); - + // Resize and initialize the // lvs lvs.reinit (n_inner, n_outer); @@ -857,7 +857,7 @@ MappingQ::apply_laplace_vector(const Table<2,double> &lvs, // reason for not aborting there // any more [WB] Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree)); - + const unsigned int n_inner_apply=lvs.n_rows(); Assert(n_inner_apply==n_inner || n_inner_apply==(degree-1)*(degree-1), ExcInternalError()); @@ -903,13 +903,13 @@ MappingQ::compute_mapping_support_points( // cell { a.resize(GeometryInfo::vertices_per_cell); - + for (unsigned int i=0; i::vertices_per_cell; ++i) a[i] = cell->vertex(i); } } - + template void MappingQ::compute_support_points_laplace(const typename Triangulation::cell_iterator &cell, @@ -920,7 +920,7 @@ MappingQ::compute_support_points_laplace(const typename Triangulat a.resize(GeometryInfo::vertices_per_cell); for (unsigned int i=0; i::vertices_per_cell; ++i) a[i] = cell->vertex(i); - + if (degree>1) switch (dim) { @@ -996,10 +996,10 @@ MappingQ<1,2>::add_line_support_points (const Triangulation<1,2>::cell_iterator // boundary description { std::vector > line_points (degree-1); - + const Boundary * const boundary = &(cell->get_triangulation().get_boundary(cell->material_id())); - + boundary->get_intermediate_points_on_line (cell, line_points); // Append all points a.insert (a.end(), line_points.begin(), line_points.end()); @@ -1028,7 +1028,7 @@ MappingQ::add_line_support_points (const typename Triangulationget_triangulation().get_boundary(cell->material_id()): &straight_boundary); - + a.push_back(boundary->get_new_point_on_line(line)); }; } @@ -1039,7 +1039,7 @@ MappingQ::add_line_support_points (const typename Triangulation > line_points (degree-1); - + // loop over each of the lines, // and if it is at the // boundary, then first get the @@ -1049,14 +1049,14 @@ MappingQ::add_line_support_points (const typename Triangulation::lines_per_cell; ++line_no) { const typename Triangulation::line_iterator line = cell->line(line_no); - + const Boundary * const boundary = (line->at_boundary()? &line->get_triangulation().get_boundary(line->boundary_indicator()) : (dim != spacedim) ? &line->get_triangulation().get_boundary(cell->material_id()) : &straight_boundary); - + boundary->get_intermediate_points_on_line (line, line_points); if (dim==3) { @@ -1073,7 +1073,7 @@ MappingQ::add_line_support_points (const typename Triangulation::cell_iterator &cell, // used if only one line of face // quad is at boundary std::vector > b(4*degree); - - + + // loop over all faces and collect // points on them for (unsigned int face_no=0; face_no::cell_iterator &cell, face_flip = cell->face_flip (face_no), face_rotation = cell->face_rotation (face_no); -#ifdef DEBUG +#ifdef DEBUG // some sanity checks up front for (unsigned int i=0; ivertex_index(i)==cell->vertex_index( @@ -1136,7 +1136,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, face_no, i, face_orientation, face_flip, face_rotation)), ExcInternalError()); #endif - + // if face at boundary, then // ask boundary object to // return intermediate points @@ -1168,7 +1168,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, for (unsigned int i=0; iline(i)->at_boundary()) ++lines_at_boundary; - + Assert(lines_at_boundary<=lines_per_face, ExcInternalError()); // if at least one of the @@ -1187,7 +1187,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // function was already // called. b.resize(4*degree); - + // b is of size // 4*degree, make sure // that this is the @@ -1195,7 +1195,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, Assert(b.size()==vertices_per_face+lines_per_face*(degree-1), ExcDimensionMismatch(b.size(), vertices_per_face+lines_per_face*(degree-1))); - + // sort the points into b. We // used access from the cell (not // from the face) to fill b, so @@ -1205,7 +1205,7 @@ add_quad_support_points(const Triangulation<3>::cell_iterator &cell, // standard orientation as well. for (unsigned int i=0; i::face_to_cell_vertices(face_no, i)]; - + for (unsigned int i=0; i::cell_iterator &cell, apply_laplace_vector(laplace_on_quad_vector, b); Assert(b.size()==4*degree+(degree-1)*(degree-1), ExcDimensionMismatch(b.size(), 4*degree+(degree-1)*(degree-1))); - + for (unsigned int i=0; i<(degree-1)*(degree-1); ++i) a.push_back(b[4*degree+i]); } @@ -1261,7 +1261,7 @@ add_quad_support_points(const Triangulation<2,3>::cell_iterator &cell, std::vector > &a) const { std::vector > quad_points ((degree-1)*(degree-1)); - + cell->get_triangulation().get_boundary(cell->material_id()) .get_intermediate_points_on_quad (cell, quad_points); for (unsigned int i=0; i::transform ( // to test further if (!q1_data->is_mapping_q1_data) { - Assert (dynamic_cast(&mapping_data) != 0, + Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast(mapping_data); // If we only use the @@ -1315,7 +1315,7 @@ MappingQ::transform ( // right tensors in it and we call // the base classes transform // function - MappingQ1::transform(input, output, *q1_data, mapping_type); + MappingQ1::transform(input, output, *q1_data, mapping_type); } @@ -1340,7 +1340,7 @@ MappingQ::transform ( // to test further if (!q1_data->is_mapping_q1_data) { - Assert (dynamic_cast(&mapping_data) != 0, + Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); const InternalData &data = static_cast(mapping_data); // If we only use the @@ -1353,7 +1353,7 @@ MappingQ::transform ( // right tensors in it and we call // the base classes transform // function - MappingQ1::transform(input, output, *q1_data, mapping_type); + MappingQ1::transform(input, output, *q1_data, mapping_type); } @@ -1373,7 +1373,7 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it std::auto_ptr mdata (dynamic_cast ( get_data(update_transformation_values, point_quadrature))); - + mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells || cell->has_boundary_lines()); @@ -1383,7 +1383,7 @@ transform_unit_to_real_cell (const typename Triangulation::cell_it &*mdata); compute_mapping_support_points(cell, p_data->mapping_support_points); - + return this->transform_unit_to_real_cell_internal(*p_data); } @@ -1398,7 +1398,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // first a Newton iteration based // on a Q1 mapping Point p_unit = MappingQ1::transform_real_to_unit_cell(cell, p); - + // then a Newton iteration based on // the full MappingQ if we need // this @@ -1412,7 +1412,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it get_data(update_transformation_values | update_transformation_gradients, point_quadrature))); - + mdata->use_mapping_q1_on_current_cell = false; std::vector > &points = mdata->mapping_support_points; @@ -1420,7 +1420,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it this->transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit); } - + return p_unit; } @@ -1442,7 +1442,7 @@ MappingQ::clone () const return new MappingQ(*this); } - + // explicit instantiation template class MappingQ; diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index 75f67acef6..ee72c8f37b 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -93,13 +93,13 @@ namespace internal void compute_shapes_virtual (const unsigned int n_shape_functions, const std::vector > &unit_points, - typename dealii::MappingQ1<1,spacedim>::InternalData& data) + typename dealii::MappingQ1<1,spacedim>::InternalData& data) { const unsigned int n_points=unit_points.size(); for (unsigned int k = 0 ; k < n_points ; ++k) { double x = unit_points[k](0); - + if (data.shape_values.size()!=0) { Assert(data.shape_values.size()==n_shape_functions*n_points, @@ -119,7 +119,7 @@ namespace internal // the following may or may not // work if dim != spacedim Assert (spacedim == 1, ExcNotImplemented()); - + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); data.second_derivative(k,0)[0][0] = 0; @@ -133,14 +133,14 @@ namespace internal void compute_shapes_virtual (const unsigned int n_shape_functions, const std::vector > &unit_points, - typename dealii::MappingQ1<2,spacedim>::InternalData& data) + typename dealii::MappingQ1<2,spacedim>::InternalData& data) { const unsigned int n_points=unit_points.size(); for (unsigned int k = 0 ; k < n_points ; ++k) { double x = unit_points[k](0); double y = unit_points[k](1); - + if (data.shape_values.size()!=0) { Assert(data.shape_values.size()==n_shape_functions*n_points, @@ -168,7 +168,7 @@ namespace internal // the following may or may not // work if dim != spacedim Assert (spacedim == 2, ExcNotImplemented()); - + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); data.second_derivative(k,0)[0][0] = 0; @@ -197,7 +197,7 @@ namespace internal void compute_shapes_virtual (const unsigned int n_shape_functions, const std::vector > &unit_points, - typename dealii::MappingQ1<3,spacedim>::InternalData& data) + typename dealii::MappingQ1<3,spacedim>::InternalData& data) { const unsigned int n_points=unit_points.size(); for (unsigned int k = 0 ; k < n_points ; ++k) @@ -205,7 +205,7 @@ namespace internal double x = unit_points[k](0); double y = unit_points[k](1); double z = unit_points[k](2); - + if (data.shape_values.size()!=0) { Assert(data.shape_values.size()==n_shape_functions*n_points, @@ -253,7 +253,7 @@ namespace internal // the following may or may not // work if dim != spacedim Assert (spacedim == 3, ExcNotImplemented()); - + Assert(data.shape_second_derivatives.size()==n_shape_functions*n_points, ExcInternalError()); data.second_derivative(k,0)[0][0] = 0; @@ -314,7 +314,7 @@ namespace internal data.second_derivative(k,5)[2][0] = (1.-y); data.second_derivative(k,6)[2][0] = -y; data.second_derivative(k,7)[2][0] = y; - + data.second_derivative(k,0)[1][2] = (1.-x); data.second_derivative(k,1)[1][2] = x; data.second_derivative(k,2)[1][2] = -(1.-x); @@ -418,7 +418,7 @@ MappingQ1::update_each (const UpdateFlags in) const if (out & (update_JxW_values | update_normal_vectors)) out |= update_boundary_forms; - + if (out & (update_covariant_transformation | update_JxW_values | update_jacobians @@ -426,7 +426,7 @@ MappingQ1::update_each (const UpdateFlags in) const | update_boundary_forms | update_normal_vectors)) out |= update_contravariant_transformation; - + if (out & (update_inverse_jacobians)) out |= update_covariant_transformation; @@ -443,7 +443,7 @@ MappingQ1::update_each (const UpdateFlags in) const out |= update_JxW_values; } - + return out; } @@ -462,7 +462,7 @@ MappingQ1::compute_data (const UpdateFlags update_flags, data.update_flags = data.update_once | data.update_each; const UpdateFlags flags(data.update_flags); - + if (flags & update_transformation_values) data.shape_values.resize(data.n_shape_functions * n_q_points); @@ -480,7 +480,7 @@ MappingQ1::compute_data (const UpdateFlags update_flags, if (flags & update_jacobian_grads) data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points); - + compute_shapes (q.get_points(), data); } @@ -512,7 +512,7 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, if (data.update_flags & update_boundary_forms) { data.aux.resize (dim-1, std::vector > (n_original_q_points)); - + // Compute tangentials to the // unit cell. const unsigned int nfaces = GeometryInfo::faces_per_cell; @@ -525,7 +525,7 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, static const int tangential_orientation[4]={-1,1,1,-1}; for (unsigned int i=0; i tang; + Tensor<1,spacedim> tang; tang[1-i/2]=tangential_orientation[i]; std::fill (data.unit_tangentials[i].begin(), data.unit_tangentials[i].end(), tang); @@ -539,7 +539,7 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, const unsigned int nd= GeometryInfo::unit_normal_direction[i]; - + // first tangential // vector in direction // of the (nd+1)%3 axis @@ -564,7 +564,7 @@ MappingQ1::compute_face_data (const UpdateFlags update_flags, } } - + template typename Mapping::InternalDataBase * @@ -607,7 +607,6 @@ MappingQ1::compute_fill (const typename Triangulation > &quadrature_points) const { - const UpdateFlags update_flags(data.current_update_flags()); // if necessary, recompute the @@ -640,7 +639,7 @@ MappingQ1::compute_fill (const typename Triangulation::compute_fill (const typename Triangulation::compute_fill (const typename Triangulation::fill_fe_values ( CellSimilarity::Similarity &cell_similarity) const { // ensure that the following cast is really correct: - Assert (dynamic_cast(&mapping_data) != 0, + Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast(mapping_data); @@ -754,7 +753,7 @@ MappingQ1::fill_fe_values ( compute_fill (cell, n_q_points, DataSetDescriptor::cell (), cell_similarity, data, quadrature_points); - + const UpdateFlags update_flags(data.current_update_flags()); const std::vector &weights=q.get_weights(); @@ -766,7 +765,7 @@ MappingQ1::fill_fe_values ( // the case <2,3> if (update_flags & (update_normal_vectors | update_JxW_values)) - { + { Assert (JxW_values.size() == n_q_points, ExcDimensionMismatch(JxW_values.size(), n_q_points)); @@ -776,7 +775,7 @@ MappingQ1::fill_fe_values ( if (cell_similarity != CellSimilarity::translation) for (unsigned int point=0; point::fill_fe_values ( data.contravariant[point]=transpose(data.contravariant[point]); JxW_values[point] = data.contravariant[point][0].norm()*weights[point]; - if(update_flags & update_normal_vectors) + if(update_flags & update_normal_vectors) { normal_vectors[point][0]=-data.contravariant[point][0][1] /data.contravariant[point][0].norm(); @@ -798,7 +797,7 @@ MappingQ1::fill_fe_values ( if ( (dim==2) && (spacedim==3) ) { data.contravariant[point]=transpose(data.contravariant[point]); - cross_product(data.contravariant[point][2], + cross_product(data.contravariant[point][2], data.contravariant[point][0], data.contravariant[point][1] ); @@ -809,17 +808,17 @@ MappingQ1::fill_fe_values ( //is stored in the 3d //subtensor of the contravariant tensor data.contravariant[point][2]/=data.contravariant[point][2].norm(); - if(update_flags & update_normal_vectors) + if(update_flags & update_normal_vectors) normal_vectors[point]=data.contravariant[point][2]; } - + } } // copy values from InternalData to vector // given by reference if (update_flags & update_jacobians) - { + { Assert (jacobians.size() == n_q_points, ExcDimensionMismatch(jacobians.size(), n_q_points)); if (cell_similarity != CellSimilarity::translation) @@ -854,7 +853,7 @@ MappingQ1::fill_fe_values ( // copy values from InternalData to vector // given by reference if (update_flags & update_inverse_jacobians) - { + { Assert (inverse_jacobians.size() == n_q_points, ExcDimensionMismatch(inverse_jacobians.size(), n_q_points)); if (cell_similarity != CellSimilarity::translation) @@ -875,10 +874,10 @@ MappingQ1<2,3>::compute_fill_face (const Triangulation<2,3>::cell_iterator &, 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()); } @@ -897,9 +896,9 @@ MappingQ1::compute_fill_face ( const DataSetDescriptor, const std::vector &, InternalData &, - std::vector > &, + std::vector > &, std::vector &, - std::vector > &, + std::vector > &, std::vector > &) const { Assert(false, ExcNotImplemented()); @@ -913,9 +912,9 @@ MappingQ1::fill_fe_face_values ( const unsigned, const Quadrature&, typename Mapping::InternalDataBase&, - std::vector >&, + std::vector >&, std::vector&, - std::vector >&, + std::vector >&, std::vector >&) const { Assert(false, ExcNotImplemented()); @@ -930,9 +929,9 @@ MappingQ1::fill_fe_subface_values ( const unsigned, const Quadrature&, typename Mapping::InternalDataBase&, - std::vector >&, + std::vector >&, std::vector&, - std::vector >&, + std::vector >&, std::vector >&) const { Assert(false, ExcNotImplemented()); @@ -952,16 +951,16 @@ MappingQ1::compute_fill_face ( const DataSetDescriptor data_set, const std::vector &weights, InternalData &data, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &boundary_forms, + std::vector > &boundary_forms, std::vector > &normal_vectors) const { compute_fill (cell, n_q_points, data_set, CellSimilarity::none, data, quadrature_points); const UpdateFlags update_flags(data.current_update_flags()); - + if (update_flags & update_boundary_forms) { Assert (boundary_forms.size()==n_q_points, @@ -971,16 +970,16 @@ MappingQ1::compute_fill_face ( ExcDimensionMismatch(normal_vectors.size(), n_q_points)); if (update_flags & update_JxW_values) Assert (JxW_values.size() == n_q_points, - ExcDimensionMismatch(JxW_values.size(), n_q_points)); - + ExcDimensionMismatch(JxW_values.size(), n_q_points)); + transform(data.unit_tangentials[face_no], data.aux[0], data, mapping_contravariant); - typename std::vector >::iterator + typename std::vector >::iterator result = boundary_forms.begin(); - const typename std::vector >::iterator + const typename std::vector >::iterator end = boundary_forms.end(); - + switch (dim) { case 2: @@ -1009,7 +1008,7 @@ MappingQ1::compute_fill_face ( default: Assert(false, ExcNotImplemented()); } - + if (update_flags & (update_normal_vectors | update_JxW_values)) for (unsigned int i=0;i::fill_fe_face_values (const typename Triangulation &q, typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &boundary_forms, + std::vector > &boundary_forms, std::vector > &normal_vectors) const { // ensure that the following cast is really correct: - Assert (dynamic_cast(&mapping_data) != 0, + Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast(mapping_data); const unsigned int n_q_points = q.size(); - + compute_fill_face (cell, face_no, deal_II_numbers::invalid_unsigned_int, n_q_points, DataSetDescriptor::face (face_no, @@ -1074,18 +1073,18 @@ MappingQ1::fill_fe_subface_values (const typename Triangulation &q, typename Mapping::InternalDataBase &mapping_data, - std::vector > &quadrature_points, + std::vector > &quadrature_points, std::vector &JxW_values, - std::vector > &boundary_forms, + std::vector > &boundary_forms, std::vector > &normal_vectors) const { // ensure that the following cast is really correct: - Assert (dynamic_cast(&mapping_data) != 0, + Assert (dynamic_cast(&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast(mapping_data); const unsigned int n_q_points = q.size(); - + compute_fill_face (cell, face_no, sub_no, n_q_points, DataSetDescriptor::subface (face_no, sub_no, @@ -1115,19 +1114,19 @@ MappingQ1::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); - Tensor<1, spacedim> auxiliary; - + Tensor<1, spacedim> auxiliary; + switch (mapping_type) { case mapping_covariant: { Assert (data.update_flags & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( } return; } - + case mapping_contravariant: { Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( } return; } - + case mapping_piola: { Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); Assert (data.update_flags & update_volume_elements, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( } return; } - + default: Assert(false, ExcNotImplemented()); } @@ -1182,20 +1181,20 @@ MappingQ1::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); Tensor<2, spacedim> aux1; Tensor<2, spacedim> aux2; - + switch (mapping_type) { case mapping_covariant: { Assert (data.update_flags & update_covariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( } return; } - + case mapping_contravariant: { Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::transform ( } return; } - + case mapping_covariant_gradient: { Assert (data.update_flags & update_contravariant_transformation, typename FEValuesBase::ExcAccessToUninitializedField()); - + for (unsigned int i=0; i::ExcAccessToUninitializedField()); 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 p_real; @@ -1365,7 +1364,7 @@ transform_real_to_unit_cell (const typename Triangulation::cell_it // perform the newton iteration. transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit); - + return p_unit; } @@ -1415,22 +1414,22 @@ transform_real_to_unit_cell_internal const unsigned int n_shapes=mdata.shape_values.size(); Assert(n_shapes!=0, ExcInternalError()); Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError()); - + std::vector > &points=mdata.mapping_support_points; Assert(points.size()==n_shapes, ExcInternalError()); - + // Newton iteration to solve // f(x)=p(x)-p=0 // x_{n+1}=x_n-[f'(x)]^{-1}f(x) - + // The start value is set to be the // center of the unit cell. - + // The shape values and derivatives // of the mapping at this point are // previously computed. - + // f(x) Point p_real(transform_unit_to_real_cell_internal(mdata)); Point f = p_real-p; @@ -1438,19 +1437,19 @@ transform_real_to_unit_cell_internal const double eps=1e-15*cell->diameter(); unsigned int loop=0; while (f.square()>eps*eps && loop++<10) - { + { // f'(x) Tensor<2,dim> df; for (unsigned int k=0; k &grad_transform=mdata.derivative(0,k); const Point &point=points[k]; - + for (unsigned int i=0; i d; Tensor<2,dim> df_1; @@ -1463,7 +1462,7 @@ transform_real_to_unit_cell_internal // shape values and derivatives // at new p_unit point compute_shapes(std::vector > (1, p_unit), mdata); - + // f(x) p_real = transform_unit_to_real_cell_internal(mdata); f = p_real-p; diff --git a/deal.II/source/hp/fe_values.cc b/deal.II/source/hp/fe_values.cc index c4b67767e9..372d4b97c0 100644 --- a/deal.II/source/hp/fe_values.cc +++ b/deal.II/source/hp/fe_values.cc @@ -309,33 +309,33 @@ namespace hp // -------------------------- FEFaceValues ------------------------- - template - FEFaceValues::FEFaceValues (const hp::MappingCollection &mapping, - const hp::FECollection &fe_collection, + template + FEFaceValues::FEFaceValues (const hp::MappingCollection &mapping, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase > (mapping, + internal::hp::FEValuesBase > (mapping, fe_collection, q_collection, update_flags) {} - template - FEFaceValues::FEFaceValues (const hp::FECollection &fe_collection, + template + FEFaceValues::FEFaceValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase > (fe_collection, + internal::hp::FEValuesBase > (fe_collection, q_collection, update_flags) {} - template + template void - FEFaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, + FEFaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index, const unsigned int mapping_index, @@ -384,9 +384,9 @@ namespace hp - template + template void - FEFaceValues::reinit (const typename dealii::DoFHandler::cell_iterator &cell, + FEFaceValues::reinit (const typename dealii::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index, const unsigned int mapping_index, @@ -425,9 +425,9 @@ namespace hp - template + template void - FEFaceValues::reinit (const typename MGDoFHandler::cell_iterator &cell, + FEFaceValues::reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index, const unsigned int mapping_index, @@ -466,9 +466,9 @@ namespace hp - template + template void - FEFaceValues::reinit (const typename Triangulation::cell_iterator &cell, + FEFaceValues::reinit (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int q_index, const unsigned int mapping_index, @@ -509,33 +509,33 @@ namespace hp // -------------------------- FESubfaceValues ------------------------- - template - FESubfaceValues::FESubfaceValues (const hp::MappingCollection &mapping, - const hp::FECollection &fe_collection, + template + FESubfaceValues::FESubfaceValues (const hp::MappingCollection &mapping, + const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase > (mapping, + internal::hp::FEValuesBase > (mapping, fe_collection, q_collection, update_flags) {} - template - FESubfaceValues::FESubfaceValues (const hp::FECollection &fe_collection, + template + FESubfaceValues::FESubfaceValues (const hp::FECollection &fe_collection, const hp::QCollection &q_collection, const UpdateFlags update_flags) : - internal::hp::FEValuesBase > (fe_collection, + internal::hp::FEValuesBase > (fe_collection, q_collection, update_flags) {} - template + template void - FESubfaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, + FESubfaceValues::reinit (const typename hp::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index, @@ -585,9 +585,9 @@ namespace hp - template + template void - FESubfaceValues::reinit (const typename dealii::DoFHandler::cell_iterator &cell, + FESubfaceValues::reinit (const typename dealii::DoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index, @@ -627,9 +627,9 @@ namespace hp - template + template void - FESubfaceValues::reinit (const typename MGDoFHandler::cell_iterator &cell, + FESubfaceValues::reinit (const typename MGDoFHandler::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index, @@ -669,9 +669,9 @@ namespace hp - template + template void - FESubfaceValues::reinit (const typename Triangulation::cell_iterator &cell, + FESubfaceValues::reinit (const typename Triangulation::cell_iterator &cell, const unsigned int face_no, const unsigned int subface_no, const unsigned int q_index, @@ -746,13 +746,12 @@ namespace internal { template class FEValuesBase >; -// not yet implemented: -// #if deal_II_dimension == 2 -// template class FEValuesBase >; -// template class FEValuesBase >; -// #endif +#if deal_II_dimension == 2 + template class FEValuesBase >; + template class FEValuesBase >; +#endif } } @@ -760,11 +759,10 @@ namespace hp { template class FEValues; -// not yet implemented: -//#if deal_II_dimension == 2 -// template class FEFaceValues; -// template class FESubfaceValues; -//#endif +#if deal_II_dimension == 2 + template class FEFaceValues; + template class FESubfaceValues; +#endif } #endif diff --git a/deal.II/source/numerics/vectors.cc b/deal.II/source/numerics/vectors.cc index c0bf0b59f7..c30a9f4c51 100644 --- a/deal.II/source/numerics/vectors.cc +++ b/deal.II/source/numerics/vectors.cc @@ -206,6 +206,101 @@ void VectorTools::interpolate_boundary_values ( ConstraintMatrix &, const std::vector &); +#if deal_II_dimension < 3 +template +void VectorTools::interpolate_boundary_values ( + const DoFHandler &, + const unsigned char, + const Function &, + std::map &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const hp::DoFHandler &, + const unsigned char, + const Function &, + std::map &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const MGDoFHandler &, + const unsigned char, + const Function &, + std::map &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const DoFHandler &, + const unsigned char, + const Function &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const hp::DoFHandler &, + const unsigned char, + const Function &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const MGDoFHandler &, + const unsigned char, + const Function &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const Mapping &, + const DoFHandler &, + const FunctionMap::type &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const Mapping &, + const hp::DoFHandler &, + const FunctionMap::type &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const Mapping &, + const MGDoFHandler &, + const FunctionMap::type &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const DoFHandler &, + const FunctionMap::type &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const hp::DoFHandler &, + const FunctionMap::type &, + ConstraintMatrix &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values ( + const MGDoFHandler &, + const FunctionMap::type &, + ConstraintMatrix &, + const std::vector &); +#endif + template void VectorTools::project_boundary_values (const Mapping &, @@ -266,22 +361,6 @@ VectorTools::compute_no_normal_flux_constraints (const DoFHandler -// (const Mapping<1> &, -// const DoFHandler<1> &, -// const unsigned char, -// const Function<1> &, -// std::map &, -// const std::vector &); -// #endif - -// the following two functions are not derived from a template in 1d -// and thus need no explicit instantiation -//#if deal_II_dimension > 1 template void VectorTools::interpolate_boundary_values (const Mapping &, @@ -299,7 +378,24 @@ void VectorTools::interpolate_boundary_values std::map &, const std::vector &); -//#endif +#if deal_II_dimension < 3 +template +void VectorTools::interpolate_boundary_values +(const Mapping &, + const DoFHandler &, + const FunctionMap::type &, + std::map &, + const std::vector &); + +template +void VectorTools::interpolate_boundary_values +(const Mapping &, + const DoFHandler &, + const unsigned char, + const Function &, + std::map &, + const std::vector &); +#endif template -- 2.39.5