<h3>deal.II</h3>
<ol>
+ <li><p>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.
+ <br>
+ (WB, 2010/11/02)
+ </p></li>
+
<li><p>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.
*
* 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 <tt>spacedim</tt> template parameter
* check the documentation of FiniteElement or the one of
* Triangulation.
class Mapping : public Subscriptor
{
public:
-
+
/**
* Virtual destructor.
*/
virtual ~Mapping ();
-
+
/**
* Transforms the point @p p on
* the unit cell to the point
transform_unit_to_real_cell (
const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<dim> &p) const = 0;
-
+
/**
* Transforms the point @p p on
* the real cell to the point
transform_real_to_unit_cell (
const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<spacedim> &p) const = 0;
-
+
/**
* Base class for internal data
* of finite element and mapping
* @p first_cell to @p true.
*/
InternalDataBase ();
-
+
/**
* Virtual destructor for
* derived classes
* by reinit.
*/
UpdateFlags update_flags;
-
+
/**
* Values computed by
* constructor.
* the first cell.
*/
virtual void clear_first_cell ();
-
+
/**
* Return an estimate (in
* bytes) or the memory
* if #update_volume_elements.
*/
std::vector<double> volume_elements;
-
+
/**
* The positions of the
* mapped (generalized)
* support points.
*/
std::vector<Point<spacedim> > support_point_values;
-
+
/*
* The Jacobian of the
* transformation in the
* points.
*/
std::vector<Tensor<2,spacedim> > support_point_gradients;
-
+
/*
* The inverse of the
* Jacobian of the
* points.
*/
std::vector<Tensor<2,spacedim> > support_point_inverse_gradients;
-
-
+
+
private:
/**
* The value returned by
*/
bool first_cell;
};
-
+
/**
* Transform a field of
* vectors accorsing to
VectorSlice<std::vector<Tensor<2,spacedim> > > output,
const InternalDataBase &internal,
const MappingType type) const = 0;
-
+
/**
* @deprecated Use transform() instead.
*/
/**
* @deprecated Use transform() instead.
*/
-
+
void
transform_contravariant (const VectorSlice<const std::vector<Tensor<2,dim> > > intput,
const unsigned int offset,
const Point<spacedim>& support_point_value(
const unsigned int index,
const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-
+
/**
* The Jacobian
* matrix of the transformation
const Tensor<2,spacedim>& support_point_gradient(
const unsigned int index,
const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-
+
/**
* The inverse Jacobian
* matrix of the transformation
const Tensor<2,spacedim>& support_point_inverse_gradient(
const unsigned int index,
const typename Mapping<dim,spacedim>::InternalDataBase &internal) const;
-
+
/**
* Return a pointer to a copy of the
* present object. The caller of this
*/
virtual
Mapping<dim,spacedim> * clone () const = 0;
-
+
/**
* Returns whether the mapping preserves
* vertex locations, i.e. whether the
DeclException0 (ExcInvalidData);
private:
-
+
/**
* Indicate fields to be updated
* in the constructor of
* See @ref UpdateFlagsEssay.
*/
virtual UpdateFlags update_once (const UpdateFlags) const = 0;
-
+
/**
* The same as update_once(),
* but for the flags to be updated for
virtual InternalDataBase*
get_face_data (const UpdateFlags flags,
const Quadrature<dim-1>& quadrature) const = 0;
-
+
/**
* Prepare internal data
* structure for transformation
const unsigned int face_no,
const Quadrature<dim-1> &quadrature,
InternalDataBase &internal,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_form,
+ std::vector<Tensor<1,spacedim> > &boundary_form,
std::vector<Point<spacedim> > &normal_vectors) const = 0;
/**
const unsigned int sub_no,
const Quadrature<dim-1> &quadrature,
InternalDataBase &internal,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_form,
+ std::vector<Tensor<1,spacedim> > &boundary_form,
std::vector<Point<spacedim> > &normal_vectors) const = 0;
/**
const typename Mapping<dim,spacedim>::InternalDataBase& internal) const
{
AssertIndexRange(index, internal.support_point_values.size());
- return internal.support_point_values[index];
+ return internal.support_point_values[index];
}
const unsigned int face_no,
const Quadrature<dim-1>& quadrature,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
+ typename std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &exterior_form,
+ typename std::vector<Tensor<1,spacedim> > &exterior_form,
typename std::vector<Point<spacedim> > &normal_vectors) const ;
/**
const unsigned int sub_no,
const Quadrature<dim-1>& quadrature,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
+ typename std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &exterior_form,
+ typename std::vector<Tensor<1,spacedim> > &exterior_form,
typename std::vector<Point<spacedim> > &normal_vectors) const ;
/**
const unsigned int face_no,
const Quadrature<dim-1> &quadrature,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
+ typename std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &boundary_form,
+ typename std::vector<Tensor<1,spacedim> > &boundary_form,
typename std::vector<Point<spacedim> > &normal_vectors) const ;
/**
const unsigned int sub_no,
const Quadrature<dim-1>& quadrature,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- typename std::vector<Point<dim> > &quadrature_points,
+ typename std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- typename std::vector<Tensor<1,dim> > &boundary_form,
+ typename std::vector<Tensor<1,spacedim> > &boundary_form,
typename std::vector<Point<spacedim> > &normal_vectors) const ;
/**
const DataSetDescriptor data_set,
const std::vector<double> &weights,
InternalData &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_form,
+ std::vector<Tensor<1,spacedim> > &boundary_form,
std::vector<Point<spacedim> > &normal_vectors) const;
/**
* @ingroup hp hpcollection
* @author Wolfgang Bangerth, 2003
*/
- template <int dim>
- class FEFaceValues : public internal::hp::FEValuesBase<dim,dim-1,dealii::FEFaceValues<dim> >
+ template <int dim, int spacedim=dim>
+ class FEFaceValues : public internal::hp::FEValuesBase<dim,dim-1,dealii::FEFaceValues<dim,spacedim> >
{
public:
/**
* <tt>DoFHandler::get_fe()</tt>
* function.
*/
- FEFaceValues (const hp::MappingCollection<dim> &mapping_collection,
- const hp::FECollection<dim> &fe_collection,
+ FEFaceValues (const hp::MappingCollection<dim,spacedim> &mapping_collection,
+ const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags);
* <tt>DoFHandler::get_fe()</tt>
* function.
*/
- FEFaceValues (const hp::FECollection<dim> &fe_collection,
+ FEFaceValues (const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags);
* this argument is specified.
*/
void
- reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
+ reinit (const typename hp::DoFHandler<dim,spacedim>::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,
* these last three arguments.
*/
void
- reinit (const typename dealii::DoFHandler<dim>::cell_iterator &cell,
+ reinit (const typename dealii::DoFHandler<dim,spacedim>::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,
* these last three arguments.
*/
void
- reinit (const typename MGDoFHandler<dim>::cell_iterator &cell,
+ reinit (const typename MGDoFHandler<dim,spacedim>::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,
* three arguments.
*/
void
- reinit (const typename Triangulation<dim>::cell_iterator &cell,
+ reinit (const typename Triangulation<dim,spacedim>::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,
* @ingroup hp hpcollection
* @author Wolfgang Bangerth, 2003
*/
- template <int dim>
- class FESubfaceValues : public internal::hp::FEValuesBase<dim,dim-1,dealii::FESubfaceValues<dim> >
+ template <int dim, int spacedim=dim>
+ class FESubfaceValues : public internal::hp::FEValuesBase<dim,dim-1,dealii::FESubfaceValues<dim,spacedim> >
{
public:
/**
* <tt>DoFHandler::get_fe()</tt>
* function.
*/
- FESubfaceValues (const hp::MappingCollection<dim> &mapping_collection,
- const hp::FECollection<dim> &fe_collection,
+ FESubfaceValues (const hp::MappingCollection<dim,spacedim> &mapping_collection,
+ const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags);
* <tt>DoFHandler::get_fe()</tt>
* function.
*/
- FESubfaceValues (const hp::FECollection<dim> &fe_collection,
+ FESubfaceValues (const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags);
* this argument is specified.
*/
void
- reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
+ reinit (const typename hp::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index = numbers::invalid_unsigned_int,
* these last three arguments.
*/
void
- reinit (const typename dealii::DoFHandler<dim>::cell_iterator &cell,
+ reinit (const typename dealii::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index = numbers::invalid_unsigned_int,
* these last three arguments.
*/
void
- reinit (const typename MGDoFHandler<dim>::cell_iterator &cell,
+ reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index = numbers::invalid_unsigned_int,
* three arguments.
*/
void
- reinit (const typename Triangulation<dim>::cell_iterator &cell,
+ reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index = numbers::invalid_unsigned_int,
const std::vector<bool> &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()),
const unsigned int n_components = DoFTools::n_components(dof);
const bool fe_is_system = (n_components != 1);
- for (typename FunctionMap<DH::space_dimension>::type::const_iterator i=function_map.begin();
+ for (typename FunctionMap<spacedim>::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));
std::vector<unsigned int> face_dofs;
face_dofs.reserve (DoFTools::max_dofs_per_face(dof));
- std::vector<Point<DH::space_dimension> > dof_locations;
+ std::vector<Point<spacedim> > dof_locations;
dof_locations.reserve (DoFTools::max_dofs_per_face(dof));
// array to store the values of
// the interpolation points of all
// finite elements that may ever be
// in use
- hp::FECollection<dim> finite_elements (dof.get_fe());
+ hp::FECollection<dim,spacedim> finite_elements (dof.get_fe());
hp::QCollection<dim-1> q_collection;
for (unsigned int f=0; f<finite_elements.size(); ++f)
{
- const FiniteElement<dim> &fe = finite_elements[f];
+ const FiniteElement<dim,spacedim> &fe = finite_elements[f];
// generate a quadrature rule
// on the face from the unit
// hp::FEFaceValues object that we
// can use to evaluate the boundary
// values at
- hp::MappingCollection<dim> mapping_collection (mapping);
- hp::FEFaceValues<dim> x_fe_values (mapping_collection, finite_elements, q_collection,
- update_quadrature_points);
+ hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
+ hp::FEFaceValues<dim,spacedim> x_fe_values (mapping_collection, finite_elements, q_collection,
+ update_quadrature_points);
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
++face_no)
{
- const FiniteElement<dim,DH::space_dimension> &fe = cell->get_fe();
+ const FiniteElement<dim,spacedim> &fe = cell->get_fe();
// we can presently deal only with
// primitive elements for boundary
{
// face is of the right component
x_fe_values.reinit(cell, face_no);
- const FEFaceValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+ const FEFaceValues<dim,spacedim> &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<Point<DH::space_dimension> > &dof_locations
+ const std::vector<Point<spacedim> > &dof_locations
= fe_values.get_quadrature_points ();
if (fe_is_system)
&& ! 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
//
FEFaceValuesBase<dim,spacedim> (quadrature.size(),
fe.dofs_per_cell,
update_flags,
- StaticMappingQ1<dim>::mapping,
+ StaticMappingQ1<dim,spacedim>::mapping,
fe, quadrature)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
FEFaceValuesBase<dim,spacedim> (quadrature.size(),
fe.dofs_per_cell,
update_flags,
- StaticMappingQ1<dim>::mapping,
+ StaticMappingQ1<dim,spacedim>::mapping,
fe, quadrature)
{
Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
template class FEValues<deal_II_dimension,deal_II_dimension+1>;
template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
CellIterator<DoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
-//template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
-// CellIterator<MGDoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
-
-// #if deal_II_dimension == 2
-// template class FEFaceValuesBase<deal_II_dimension,deal_II_dimension+1>;
-// template class FEFaceValues<deal_II_dimension,deal_II_dimension+1>;
-// template class FESubfaceValues<deal_II_dimension,deal_II_dimension+1>;
-// #endif
+template class FEValuesBase<deal_II_dimension,deal_II_dimension+1>::
+ CellIterator<MGDoFHandler<deal_II_dimension,deal_II_dimension+1>::cell_iterator>;
+
+#if deal_II_dimension == 2
+template class FEFaceValuesBase<deal_II_dimension,deal_II_dimension+1>;
+template class FEFaceValues<deal_II_dimension,deal_II_dimension+1>;
+template class FESubfaceValues<deal_II_dimension,deal_II_dimension+1>;
+#endif
#endif
template<int dim, int spacedim>
unsigned int
-MappingQ<dim,spacedim>::InternalData::memory_consumption () const
+MappingQ<dim,spacedim>::InternalData::memory_consumption () const
{
return (MappingQ1<dim,spacedim>::InternalData::memory_consumption () +
MemoryConsumption::memory_consumption (unit_normals) +
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)
{
ExcInternalError());
grads.resize(n_shape_functions);
}
-
+
// // dummy variable of size 0
std::vector<Tensor<2,dim> > grad2;
if (data.shape_second_derivatives.size()!=0)
grad2.resize(n_shape_functions);
}
-
+
if (data.shape_values.size()!=0 || data.shape_derivatives.size()!=0)
for (unsigned int point=0; point<n_points; ++point)
{
tensor_pols->compute(unit_points[point], values, grads, grad2);
-
+
if (data.shape_values.size()!=0)
for (unsigned int i=0; i<n_shape_functions; ++i)
data.shape(point,renumber[i]) = values[i];
-
+
if (data.shape_derivatives.size()!=0)
for (unsigned int i=0; i<n_shape_functions; ++i)
data.derivative(point,renumber[i]) = grads[i];
if (get_degree() > 1)
cell_similarity = CellSimilarity::invalid_next_cell;
}
-
+
MappingQ1<dim,spacedim>::fill_fe_values(cell, q, *p_data,
quadrature_points, JxW_values,
jacobians, jacobian_grads, inverse_jacobians,
const unsigned int face_no,
const Quadrature<dim-1> &q,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &exterior_forms,
+ std::vector<Tensor<1,spacedim> > &exterior_forms,
std::vector<Point<spacedim> > &normal_vectors) const
{
// convert data object to internal
Assert (dynamic_cast<InternalData*> (&mapping_data) != 0,
ExcInternalError());
InternalData &data = static_cast<InternalData&> (mapping_data);
-
+
// check whether this cell needs
// the full mapping or can be
// treated by a reduced Q1 mapping,
const unsigned int sub_no,
const Quadrature<dim-1> &q,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &exterior_forms,
+ std::vector<Tensor<1,spacedim> > &exterior_forms,
std::vector<Point<spacedim> > &normal_vectors) const
{
// convert data object to internal
// 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
Assert (sizeof(loqv2)/sizeof(loqv2[0]) ==
n_inner_2d * n_outer_2d,
ExcInternalError());
-
+
break;
}
// (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.,
Assert (sizeof(loqv3)/sizeof(loqv3[0]) ==
n_inner_2d * n_outer_2d,
ExcInternalError());
-
+
loqv_ptr=&loqv3[0];
-
+
break;
}
0.2231273865431891, 0.1346851306015187,
0.03812914216116723, 0.02913160002633253,
0.02200737428129391, 0.01600835564431222,
-
+
0.00664803151334206, 0.006648031513342719,
0.002873452861657458, 0.002873452861657626,
0.07903572682584378, 0.05969238281250031,
0.03812914216116729, 0.1346851306015185,
0.2231273865431898, 0.01600835564431217,
0.02200737428129394, 0.02913160002633262,
-
+
0.006648031513342073, 0.002873452861657473,
0.006648031513342726, 0.002873452861657636,
0.1527716818820238, 0.2348152760709273,
0.07903572682584376, 0.05969238281250026,
0.03619864817415824, 0.07903572682584187,
0.0596923828124998, 0.0361986481741581,
-
+
0.01106770833333302, 0.01106770833333336,
0.01106770833333337, 0.01106770833333374,
0.06770833333333424, 0.1035156250000011,
0.06770833333333422, 0.1035156250000009,
0.06770833333333436, 0.0677083333333337,
0.1035156249999988, 0.0677083333333339,
-
+
0.002873452861657185, 0.006648031513342362,
0.002873452861657334, 0.006648031513343038,
0.02496269311797779, 0.04081948955407401,
0.03619864817415819, 0.05969238281250028,
0.07903572682584407, 0.03619864817415804,
0.05969238281249986, 0.0790357268258422,
-
+
-0.001075744628906913, 0.00191429223907134,
0.07405921850311592, -0.001075744628905865,
0.03812914216116729, 0.1346851306015185,
0.02496269311797776, 0.04081948955407392,
0.02496269311797785, 0.1527716818820233,
0.2348152760709258, 0.1527716818820236,
-
+
0.001914292239071237, -0.001075744628906803,
-0.001075744628906778, 0.07405921850311617,
0.01600835564431228, 0.02200737428129401,
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
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
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
// points on the unit cell
const QGauss<dim> 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<double> S(n_inner);
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.
T(i,k) += contract(quadrature_data.derivative(point, n_outer+i),
quadrature_data.derivative(point, k))
*quadrature.weight(point);
-
+
FullMatrix<double> S_1(n_inner);
S_1.invert(S);
-
+
FullMatrix<double> 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);
// 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());
// cell
{
a.resize(GeometryInfo<dim>::vertices_per_cell);
-
+
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
a[i] = cell->vertex(i);
}
}
-
+
template<int dim, int spacedim>
void
MappingQ<dim,spacedim>::compute_support_points_laplace(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
a.resize(GeometryInfo<dim>::vertices_per_cell);
for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
a[i] = cell->vertex(i);
-
+
if (degree>1)
switch (dim)
{
// boundary description
{
std::vector<Point<spacedim> > line_points (degree-1);
-
+
const Boundary<dim,spacedim> * 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());
(dim != spacedim) ?
&line->get_triangulation().get_boundary(cell->material_id()):
&straight_boundary);
-
+
a.push_back(boundary->get_new_point_on_line(line));
};
}
// boundary description
{
std::vector<Point<spacedim> > line_points (degree-1);
-
+
// loop over each of the lines,
// and if it is at the
// boundary, then first get the
for (unsigned int line_no=0; line_no<GeometryInfo<dim>::lines_per_cell; ++line_no)
{
const typename Triangulation<dim,spacedim>::line_iterator line = cell->line(line_no);
-
+
const Boundary<dim,spacedim> * 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)
{
// correct orientation. simply
// append all points
a.insert (a.end(), line_points.begin(), line_points.end());
-
+
}
}
}
// used if only one line of face
// quad is at boundary
std::vector<Point<3> > b(4*degree);
-
-
+
+
// loop over all faces and collect
// points on them
for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
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; i<vertices_per_face; ++i)
Assert(face->vertex_index(i)==cell->vertex_index(
face_no, i, face_orientation, face_flip, face_rotation)),
ExcInternalError());
#endif
-
+
// if face at boundary, then
// ask boundary object to
// return intermediate points
for (unsigned int i=0; i<lines_per_face; ++i)
if (face->line(i)->at_boundary())
++lines_at_boundary;
-
+
Assert(lines_at_boundary<=lines_per_face, ExcInternalError());
// if at least one of the
// function was already
// called.
b.resize(4*degree);
-
+
// b is of size
// 4*degree, make sure
// that this is the
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
// standard orientation as well.
for (unsigned int i=0; i<vertices_per_face; ++i)
b[i]=a[GeometryInfo<3>::face_to_cell_vertices(face_no, i)];
-
+
for (unsigned int i=0; i<lines_per_face; ++i)
for (unsigned int j=0; j<degree-1; ++j)
b[vertices_per_face+i*(degree-1)+j]=
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]);
}
std::vector<Point<3> > &a) const
{
std::vector<Point<3> > 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<quad_points.size(); ++i)
// to test further
if (!q1_data->is_mapping_q1_data)
{
- Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData&>(mapping_data);
// If we only use the
// right tensors in it and we call
// the base classes transform
// function
- MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
+ MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
}
// to test further
if (!q1_data->is_mapping_q1_data)
{
- Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData&>(mapping_data);
// If we only use the
// right tensors in it and we call
// the base classes transform
// function
- MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
+ MappingQ1<dim,spacedim>::transform(input, output, *q1_data, mapping_type);
}
std::auto_ptr<InternalData>
mdata (dynamic_cast<InternalData *> (
get_data(update_transformation_values, point_quadrature)));
-
+
mdata->use_mapping_q1_on_current_cell = !(use_mapping_q_on_all_cells ||
cell->has_boundary_lines());
&*mdata);
compute_mapping_support_points(cell, p_data->mapping_support_points);
-
+
return this->transform_unit_to_real_cell_internal(*p_data);
}
// first a Newton iteration based
// on a Q1 mapping
Point<dim> p_unit = MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
-
+
// then a Newton iteration based on
// the full MappingQ if we need
// this
get_data(update_transformation_values |
update_transformation_gradients,
point_quadrature)));
-
+
mdata->use_mapping_q1_on_current_cell = false;
std::vector<Point<spacedim> > &points = mdata->mapping_support_points;
this->transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
}
-
+
return p_unit;
}
return new MappingQ<dim,spacedim>(*this);
}
-
+
// explicit instantiation
template class MappingQ<deal_II_dimension>;
void
compute_shapes_virtual (const unsigned int n_shape_functions,
const std::vector<Point<1> > &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,
// 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;
void
compute_shapes_virtual (const unsigned int n_shape_functions,
const std::vector<Point<2> > &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,
// 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;
void
compute_shapes_virtual (const unsigned int n_shape_functions,
const std::vector<Point<3> > &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)
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,
// 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;
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);
if (out & (update_JxW_values
| update_normal_vectors))
out |= update_boundary_forms;
-
+
if (out & (update_covariant_transformation
| update_JxW_values
| update_jacobians
| update_boundary_forms
| update_normal_vectors))
out |= update_contravariant_transformation;
-
+
if (out & (update_inverse_jacobians))
out |= update_covariant_transformation;
out |= update_JxW_values;
}
-
+
return out;
}
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);
if (flags & update_jacobian_grads)
data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
-
+
compute_shapes (q.get_points(), data);
}
if (data.update_flags & update_boundary_forms)
{
data.aux.resize (dim-1, std::vector<Tensor<1,spacedim> > (n_original_q_points));
-
+
// Compute tangentials to the
// unit cell.
const unsigned int nfaces = GeometryInfo<dim>::faces_per_cell;
static const int tangential_orientation[4]={-1,1,1,-1};
for (unsigned int i=0; i<nfaces; ++i)
{
- Tensor<1,spacedim> 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);
const unsigned int nd=
GeometryInfo<dim>::unit_normal_direction[i];
-
+
// first tangential
// vector in direction
// of the (nd+1)%3 axis
}
}
-
+
template<int dim, int spacedim>
typename Mapping<dim,spacedim>::InternalDataBase *
InternalData &data,
std::vector<Point<spacedim> > &quadrature_points) const
{
-
const UpdateFlags update_flags(data.current_update_flags());
// if necessary, recompute the
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int k=0; k<data.n_shape_functions; ++k)
- quadrature_points[point]
+ quadrature_points[point]
+= data.shape(point+data_set,k) * data.mapping_support_points[k];
}
else if (dim == spacedim - 1)
{
- // CODIMENSION 1
+ // CODIMENSION 1
// calculate left-inversion of the
// contravariant matrices to obtain
// covariant ones (auxiliary
for (unsigned int point=0; point<n_q_points; ++point)
{
contravariant_matrix.
- copy_from(data.contravariant[point],0,spacedim-1,0,dim-1);
+ copy_from(data.contravariant[point],0,spacedim-1,0,dim-1);
covariant_matrix.
left_invert(contravariant_matrix);
covariant_matrix.
CellSimilarity::Similarity &cell_similarity) const
{
// ensure that the following cast is really correct:
- Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
ExcInternalError());
InternalData &data = static_cast<InternalData&>(mapping_data);
compute_fill (cell, n_q_points, DataSetDescriptor::cell (), cell_similarity,
data, quadrature_points);
-
+
const UpdateFlags update_flags(data.current_update_flags());
const std::vector<double> &weights=q.get_weights();
// 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));
if (cell_similarity != CellSimilarity::translation)
for (unsigned int point=0; point<n_q_points; ++point)
- {
+ {
if (dim==spacedim)
JxW_values[point]
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();
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]
);
//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)
// 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)
const DataSetDescriptor,
const std::vector<double>&,
InternalData &,
- std::vector<Point<2> >&,
+ std::vector<Point<3> >&,
std::vector<double>&,
- std::vector<Tensor<1,2> > &,
- std::vector<Point<3> > &) const
+ std::vector<Tensor<1,3> > &,
+ std::vector<Point<3> > &) const
{
Assert(false, ExcNotImplemented());
}
const DataSetDescriptor,
const std::vector<double> &,
InternalData &,
- std::vector<Point<dim> > &,
+ std::vector<Point<spacedim> > &,
std::vector<double> &,
- std::vector<Tensor<1,dim> > &,
+ std::vector<Tensor<1,spacedim> > &,
std::vector<Point<spacedim> > &) const
{
Assert(false, ExcNotImplemented());
const unsigned,
const Quadrature<dim-1>&,
typename Mapping<dim,spacedim>::InternalDataBase&,
- std::vector<Point<dim> >&,
+ std::vector<Point<spacedim> >&,
std::vector<double>&,
- std::vector<Tensor<1,dim> >&,
+ std::vector<Tensor<1,spacedim> >&,
std::vector<Point<spacedim> >&) const
{
Assert(false, ExcNotImplemented());
const unsigned,
const Quadrature<dim-1>&,
typename Mapping<dim,spacedim>::InternalDataBase&,
- std::vector<Point<dim> >&,
+ std::vector<Point<spacedim> >&,
std::vector<double>&,
- std::vector<Tensor<1,dim> >&,
+ std::vector<Tensor<1,spacedim> >&,
std::vector<Point<spacedim> >&) const
{
Assert(false, ExcNotImplemented());
const DataSetDescriptor data_set,
const std::vector<double> &weights,
InternalData &data,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_forms,
+ std::vector<Tensor<1,spacedim> > &boundary_forms,
std::vector<Point<spacedim> > &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,
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<Tensor<1,dim> >::iterator
+ typename std::vector<Tensor<1,spacedim> >::iterator
result = boundary_forms.begin();
- const typename std::vector<Tensor<1,dim> >::iterator
+ const typename std::vector<Tensor<1,spacedim> >::iterator
end = boundary_forms.end();
-
+
switch (dim)
{
case 2:
default:
Assert(false, ExcNotImplemented());
}
-
+
if (update_flags & (update_normal_vectors
| update_JxW_values))
for (unsigned int i=0;i<boundary_forms.size();++i)
const unsigned int face_no,
const Quadrature<dim-1> &q,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_forms,
+ std::vector<Tensor<1,spacedim> > &boundary_forms,
std::vector<Point<spacedim> > &normal_vectors) const
{
// ensure that the following cast is really correct:
- Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
ExcInternalError());
InternalData &data = static_cast<InternalData&>(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,
const unsigned int sub_no,
const Quadrature<dim-1> &q,
typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<spacedim> > &quadrature_points,
std::vector<double> &JxW_values,
- std::vector<Tensor<1,dim> > &boundary_forms,
+ std::vector<Tensor<1,spacedim> > &boundary_forms,
std::vector<Point<spacedim> > &normal_vectors) const
{
// ensure that the following cast is really correct:
- Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
ExcInternalError());
InternalData &data = static_cast<InternalData&>(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,
const MappingType mapping_type) const
{
AssertDimension (input.size(), output.size());
- Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData&>(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<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d=0;d<dim;++d)
}
return;
}
-
+
case mapping_contravariant:
{
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d=0;d<dim;++d)
}
return;
}
-
+
case mapping_piola:
{
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
Assert (data.update_flags & update_volume_elements,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d=0;d<dim;++d)
}
return;
}
-
+
default:
Assert(false, ExcNotImplemented());
}
const MappingType mapping_type) const
{
AssertDimension (input.size(), output.size());
- Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+ Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
ExcInternalError());
const InternalData &data = static_cast<const InternalData&>(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<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d1=0;d1<dim;++d1)
}
return;
}
-
+
case mapping_contravariant:
{
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d1=0;d1<dim;++d1)
}
return;
}
-
+
case mapping_covariant_gradient:
{
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d1=0;d1<dim;++d1)
for (unsigned int d2=0;d2<dim;++d2)
aux1[d1][d2] = input[i][d1][d2];
-
+
contract(aux2, aux1, data.covariant[i]);
contract(output[i], data.covariant[i], aux2);
}
-
+
return;
}
-
+
case mapping_contravariant_gradient:
{
Assert (data.update_flags & update_covariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d1=0;d1<dim;++d1)
for (unsigned int d2=0;d2<dim;++d2)
aux1[d1][d2] = input[i][d1][d2];
-
+
contract(aux2, aux1, data.covariant[i]);
contract(output[i], data.contravariant[i], aux2);
}
-
+
return;
}
-
+
case mapping_piola:
{
Assert (data.update_flags & update_covariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
Assert (data.update_flags & update_volume_elements,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
-
+
for (unsigned int i=0; i<output.size(); ++i)
{
for (unsigned int d1=0;d1<dim;++d1)
for (unsigned int d2=0;d2<dim;++d2)
aux1[d1][d2] = input[i][d1][d2] / data.volume_elements[i];
-
+
contract (aux2, aux1, data.covariant[i]);
contract (output[i], data.contravariant[i], aux2);
}
{
const unsigned int n_mapping_points=data.mapping_support_points.size();
Assert(data.shape_values.size()==n_mapping_points, ExcInternalError());
-
+
// use now the InternalData to
// compute the point in real space.
Point<spacedim> p_real;
// perform the newton iteration.
transform_real_to_unit_cell_internal(cell, p, *mdata, p_unit);
-
+
return p_unit;
}
const unsigned int n_shapes=mdata.shape_values.size();
Assert(n_shapes!=0, ExcInternalError());
Assert(mdata.shape_derivatives.size()==n_shapes, ExcInternalError());
-
+
std::vector<Point<spacedim> > &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<spacedim> p_real(transform_unit_to_real_cell_internal(mdata));
Point<spacedim> f = p_real-p;
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<mdata.n_shape_functions; ++k)
{
const Tensor<1,dim> &grad_transform=mdata.derivative(0,k);
const Point<spacedim> &point=points[k];
-
+
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
df[i][j]+=point[i]*grad_transform[j];
}
-
+
// Solve [f'(x)]d=f(x)
Tensor<1,dim> d;
Tensor<2,dim> df_1;
// shape values and derivatives
// at new p_unit point
compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
-
+
// f(x)
p_real = transform_unit_to_real_cell_internal(mdata);
f = p_real-p;
// -------------------------- FEFaceValues -------------------------
- template <int dim>
- FEFaceValues<dim>::FEFaceValues (const hp::MappingCollection<dim> &mapping,
- const hp::FECollection<dim> &fe_collection,
+ template <int dim, int spacedim>
+ FEFaceValues<dim,spacedim>::FEFaceValues (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags)
:
- internal::hp::FEValuesBase<dim,dim-1,dealii::FEFaceValues<dim> > (mapping,
+ internal::hp::FEValuesBase<dim,dim-1,dealii::FEFaceValues<dim,spacedim> > (mapping,
fe_collection,
q_collection,
update_flags)
{}
- template <int dim>
- FEFaceValues<dim>::FEFaceValues (const hp::FECollection<dim> &fe_collection,
+ template <int dim, int spacedim>
+ FEFaceValues<dim,spacedim>::FEFaceValues (const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags)
:
- internal::hp::FEValuesBase<dim,dim-1,dealii::FEFaceValues<dim> > (fe_collection,
+ internal::hp::FEValuesBase<dim,dim-1,dealii::FEFaceValues<dim,spacedim> > (fe_collection,
q_collection,
update_flags)
{}
- template <int dim>
+ template <int dim, int spacedim>
void
- FEFaceValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
+ FEFaceValues<dim,spacedim>::reinit (const typename hp::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int q_index,
const unsigned int mapping_index,
- template <int dim>
+ template <int dim, int spacedim>
void
- FEFaceValues<dim>::reinit (const typename dealii::DoFHandler<dim>::cell_iterator &cell,
+ FEFaceValues<dim,spacedim>::reinit (const typename dealii::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int q_index,
const unsigned int mapping_index,
- template <int dim>
+ template <int dim, int spacedim>
void
- FEFaceValues<dim>::reinit (const typename MGDoFHandler<dim>::cell_iterator &cell,
+ FEFaceValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int q_index,
const unsigned int mapping_index,
- template <int dim>
+ template <int dim, int spacedim>
void
- FEFaceValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
+ FEFaceValues<dim,spacedim>::reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int q_index,
const unsigned int mapping_index,
// -------------------------- FESubfaceValues -------------------------
- template <int dim>
- FESubfaceValues<dim>::FESubfaceValues (const hp::MappingCollection<dim> &mapping,
- const hp::FECollection<dim> &fe_collection,
+ template <int dim, int spacedim>
+ FESubfaceValues<dim,spacedim>::FESubfaceValues (const hp::MappingCollection<dim,spacedim> &mapping,
+ const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags)
:
- internal::hp::FEValuesBase<dim,dim-1,dealii::FESubfaceValues<dim> > (mapping,
+ internal::hp::FEValuesBase<dim,dim-1,dealii::FESubfaceValues<dim,spacedim> > (mapping,
fe_collection,
q_collection,
update_flags)
{}
- template <int dim>
- FESubfaceValues<dim>::FESubfaceValues (const hp::FECollection<dim> &fe_collection,
+ template <int dim, int spacedim>
+ FESubfaceValues<dim,spacedim>::FESubfaceValues (const hp::FECollection<dim,spacedim> &fe_collection,
const hp::QCollection<dim-1> &q_collection,
const UpdateFlags update_flags)
:
- internal::hp::FEValuesBase<dim,dim-1,dealii::FESubfaceValues<dim> > (fe_collection,
+ internal::hp::FEValuesBase<dim,dim-1,dealii::FESubfaceValues<dim,spacedim> > (fe_collection,
q_collection,
update_flags)
{}
- template <int dim>
+ template <int dim, int spacedim>
void
- FESubfaceValues<dim>::reinit (const typename hp::DoFHandler<dim>::cell_iterator &cell,
+ FESubfaceValues<dim,spacedim>::reinit (const typename hp::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index,
- template <int dim>
+ template <int dim, int spacedim>
void
- FESubfaceValues<dim>::reinit (const typename dealii::DoFHandler<dim>::cell_iterator &cell,
+ FESubfaceValues<dim,spacedim>::reinit (const typename dealii::DoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index,
- template <int dim>
+ template <int dim, int spacedim>
void
- FESubfaceValues<dim>::reinit (const typename MGDoFHandler<dim>::cell_iterator &cell,
+ FESubfaceValues<dim,spacedim>::reinit (const typename MGDoFHandler<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index,
- template <int dim>
+ template <int dim, int spacedim>
void
- FESubfaceValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
+ FESubfaceValues<dim,spacedim>::reinit (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const unsigned int q_index,
{
template class FEValuesBase<deal_II_dimension,deal_II_dimension,
dealii::FEValues<deal_II_dimension,deal_II_dimension+1> >;
-// not yet implemented:
-// #if deal_II_dimension == 2
-// template class FEValuesBase<deal_II_dimension,deal_II_dimension-1,
-// dealii::FEFaceValues<deal_II_dimension> >;
-// template class FEValuesBase<deal_II_dimension,deal_II_dimension-1,
-// dealii::FESubfaceValues<deal_II_dimension> >;
-// #endif
+#if deal_II_dimension == 2
+ template class FEValuesBase<deal_II_dimension,deal_II_dimension-1,
+ dealii::FEFaceValues<deal_II_dimension,deal_II_dimension+1> >;
+ template class FEValuesBase<deal_II_dimension,deal_II_dimension-1,
+ dealii::FESubfaceValues<deal_II_dimension,deal_II_dimension+1> >;
+#endif
}
}
{
template class FEValues<deal_II_dimension, deal_II_dimension+1>;
-// not yet implemented:
-//#if deal_II_dimension == 2
-// template class FEFaceValues<deal_II_dimension, deal_II_dimension+1>;
-// template class FESubfaceValues<deal_II_dimension, deal_II_dimension+1>;
-//#endif
+#if deal_II_dimension == 2
+ template class FEFaceValues<deal_II_dimension, deal_II_dimension+1>;
+ template class FESubfaceValues<deal_II_dimension, deal_II_dimension+1>;
+#endif
}
#endif
ConstraintMatrix &,
const std::vector<bool> &);
+#if deal_II_dimension < 3
+template
+void VectorTools::interpolate_boundary_values (
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ std::map<unsigned int,double> &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ std::map<unsigned int,double> &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const MGDoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ std::map<unsigned int,double> &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const MGDoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const Mapping<deal_II_dimension,deal_II_dimension+1> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const Mapping<deal_II_dimension,deal_II_dimension+1> &,
+ const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const Mapping<deal_II_dimension,deal_II_dimension+1> &,
+ const MGDoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const hp::DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values (
+ const MGDoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ ConstraintMatrix &,
+ const std::vector<bool> &);
+#endif
+
template
void VectorTools::project_boundary_values<deal_II_dimension>
(const Mapping<deal_II_dimension> &,
-// // Due to introducing the DoFHandler as a template parameter,
-// // the following instantiations are required in 1d
-// #if deal_II_dimension == 1
-// template
-// void VectorTools::interpolate_boundary_values<deal_II_dimension>
-// (const Mapping<1> &,
-// const DoFHandler<1> &,
-// const unsigned char,
-// const Function<1> &,
-// std::map<unsigned int,double> &,
-// const std::vector<bool> &);
-// #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<deal_II_dimension> &,
std::map<unsigned int,double> &,
const std::vector<bool> &);
-//#endif
+#if deal_II_dimension < 3
+template
+void VectorTools::interpolate_boundary_values
+(const Mapping<deal_II_dimension,deal_II_dimension+1> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const FunctionMap<deal_II_dimension+1>::type &,
+ std::map<unsigned int,double> &,
+ const std::vector<bool> &);
+
+template
+void VectorTools::interpolate_boundary_values
+(const Mapping<deal_II_dimension,deal_II_dimension+1> &,
+ const DoFHandler<deal_II_dimension,deal_II_dimension+1> &,
+ const unsigned char,
+ const Function<deal_II_dimension+1> &,
+ std::map<unsigned int,double> &,
+ const std::vector<bool> &);
+#endif
template