// $Id$
// Version: $Name$
//
-// Copyright (C) 2005, 2006 by the deal.II authors
+// Copyright (C) 2005, 2006, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* in #mapping_type. Therefore, each constructor should contain a line
* like:
* @verbatim
- * this->mapping_type = this->independent_on_cartesian;
+ * this->mapping_type = this->mapping_none;
* @endverbatim
*
* @see PolynomialsBDM, PolynomialsRaviartThomas
virtual UpdateFlags update_each (const UpdateFlags flags) const;
protected:
- /**
- * Different options for
- * transforming the basis
- * functions from the reference
- * cell to the actual mesh cell.
- *
- * Most vector valued elements
- * either transform shape
- * functions to keep node values
- * on edges meaningful. Still, in
- * special cases, it may be
- * possible to avoid the mapping.
- */
- enum MappingType {
- /**
- * No mapping has been
- * selected, throw an error
- * if needed.
- */
- no_mapping,
- /**
- * Shape functions do not
- * depend on actual mesh
- * cell
- */
- independent,
- /**
- * Shape functions do not
- * depend on actual mesh
- * cell. The mapping class
- * must be
- * MappingCartesian.
- */
- independent_on_cartesian,
- /**
- * Shape functions are
- * transformed covariant.
- */
- covariant,
- /**
- * Shape functions are
- * transformed
- * contravariant.
- */
- contravariant
- };
-
/**
* The mapping type to be used to
* map shape functions from the
- * reference cell to the msh
+ * reference cell to the mesh
* cell.
*/
MappingType mapping_type;
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* the Mapping.
*/
update_transformation_gradients = 0x2000,
-
+ //! Determinant of the Jacobian
+ /**
+ * Compute the volume element
+ * in each quadrature point.
+ */
+ update_volume_elements = 0x4000,
/**
* Compute the JxW values
* on faces for the cell mapping
* is used in conjunction with
* H_div subspaces like RT and ABF.
*/
- update_cell_JxW_values = 0x4000,
+ update_cell_JxW_values = 0x8000,
/**
* Update the location of the
* mapped generalized support
* points of the element.
*/
- update_support_points = 0x8000,
+ update_support_points = 0x10000,
/**
* Update the Jacobian of the
* mapping in generalized
* support points.
*/
- update_support_jacobians = 0x10000,
+ update_support_jacobians = 0x20000,
/**
* Update the inverse Jacobian
* of the mapping in
* generalized support points.
*/
- update_support_inverse_jacobians = 0x20000,
+ update_support_inverse_jacobians = 0x40000,
/**
* @deprecated Update
* quadrature points
* vectors. Use
* update_face_normal_vectors
*/
- update_normal_vectors = update_face_normal_vectors
-
+ update_normal_vectors = update_face_normal_vectors,
+ //! Values needed for Piola transform
+ /**
+ * Combination of the flags
+ * needed for Piola transform
+ * of Hdiv elements.
+ */
+ update_piola = update_volume_elements | update_contravariant_transformation
};
template <int dim, int spacedim> class FEFaceValues;
template <int dim, int spacedim> class FESubfaceValues;
-//TODO: Offset in transform functions should be replaced by initializing VectorSlice correctly
-
/**
* The transformation type used
* for the Mapping::transform() functions.
* for vectors. If such a
* MappingType is applied to a
* rank 2 tensor, it is implied
- * that the resulting Tensor
- * corresponds to the derivative
- * of the vector.
+ * that the mapping is applied to
+ * each column.
*/
enum MappingType
{
+/// No mapping
+ mapping_none = 0x0000,
/// Covariant mapping
mapping_covariant = 0x0001,
/// Contravariant mapping
mapping_contravariant = 0x0002,
/// The Piola transform usually used for Hdiv elements
- mapping_piola = 0x0003,
+ mapping_piola = 0x0100,
+/// The Piola transform for the derivative of an Hdiv element
+ mapping_piola_gradient = 0x0101,
/// The mapping used for Raviart-Thomas elements
mapping_raviart_thomas = mapping_piola,
/// The mapping used for BDM elements
* object.
*/
virtual unsigned int memory_consumption () const;
+
+ /**
+ * The determinant of the
+ * Jacobian in each
+ * quadrature point. Filled
+ * if #update_volume_elements.
+ */
+ std::vector<double> volume_elements;
/**
* The positions of the
*/
Tensor<1,dim> length;
+ /**
+ * The volume element
+ */
+ double volume_element;
+
/**
* Vector of all quadrature
* points. Especially, all
// $Id$
// Version: $Name$
//
-// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Assert (dim >= 2, ExcImpossibleInDim(dim));
const unsigned int n_dofs = this->dofs_per_cell;
- this->mapping_type = this->contravariant;
+ this->mapping_type = mapping_piola;
// First, initialize the
// generalized support points and
// quadrature weights, since they
if (flags & update_values)
out |= update_values | update_covariant_transformation
- | update_contravariant_transformation
+ | update_piola
| update_cell_JxW_values
| update_JxW_values;
if (flags & update_gradients)
out |= update_gradients | update_covariant_transformation
- | update_contravariant_transformation
+ | update_piola
| update_cell_JxW_values
| update_JxW_values;
//TODO: Set update flags appropriately and figure out, how the second
template <class POLY, int dim, int spacedim>
typename Mapping<dim,spacedim>::InternalDataBase *
-FE_PolyTensor<POLY,dim,spacedim>::get_data (const UpdateFlags update_flags,
- const Mapping<dim,spacedim> &mapping,
- const Quadrature<dim> &quadrature) const
+FE_PolyTensor<POLY,dim,spacedim>::get_data (
+ const UpdateFlags update_flags,
+ const Mapping<dim,spacedim> &mapping,
+ const Quadrature<dim> &quadrature) const
{
// generate a new data object and
// initialize some fields
// possible
InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
-// Assert(mapping_type == independent
-// || ( mapping_type == independent_on_cartesian
-// && dynamic_cast<const MappingCartesian<dim>*>(&mapping) != 0),
-// ExcNotImplemented());
-
const unsigned int n_q_points = quadrature.size();
const UpdateFlags flags(fe_data.current_update_flags());
{
switch (mapping_type)
{
- case independent:
- case independent_on_cartesian:
+ case mapping_none:
{
for (unsigned int k=0; k<n_q_points; ++k)
for (unsigned int d=0; d<dim; ++d)
break;
}
- case covariant:
- case contravariant:
+ case mapping_covariant:
+ case mapping_contravariant:
{
// Use auxiliary vector for
// transformation
std::vector<Tensor<1,dim> > shape_values (n_q_points);
- if (mapping_type == covariant)
+ if (mapping_type == mapping_covariant)
mapping.transform(fe_data.shape_values[i],
shape_values, mapping_data, mapping_covariant);
else
// then copy over to target:
for (unsigned int k=0; k<n_q_points; ++k)
{
+//TODO: [GK] Remove this after change to Piola complete
+
// Recompute determinant
+
const double
- J = (mapping_type == contravariant ?
+ J = (mapping_type == mapping_contravariant ?
data.JxW_values[k] / quadrature.weight(k)
:
1.0);
}
break;
}
+ case mapping_piola:
+ if (true)
+ {
+ std::vector<Tensor<1,dim> > shape_values (n_q_points);
+ mapping.transform(fe_data.shape_values[i], shape_values,
+ mapping_data, mapping_piola);
+ for (unsigned int k=0; k<n_q_points; ++k)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_values(first+d,k)
+ = sign_change[i] * shape_values[k][d];
+ break;
+ }
default:
Assert(false, ExcNotImplemented());
switch (mapping_type)
{
- case independent:
- case independent_on_cartesian:
+ case mapping_none:
{
mapping.transform(fe_data.shape_grads[i], shape_grads1,
mapping_data, mapping_covariant);
break;
}
- case covariant:
+ case mapping_covariant:
{
mapping.transform(fe_data.shape_grads[i], shape_grads1,
mapping_data, mapping_covariant);
break;
}
-
- case contravariant:
+ case mapping_piola:
+ case mapping_contravariant:
{
mapping.transform(fe_data.shape_grads[i], shape_grads1,
mapping_data, mapping_covariant);
for (unsigned int d=0; d<dim; ++d)
{
// Recompute determinant
- const double
- J = (mapping_type == contravariant ?
- data.JxW_values[k] / quadrature.weight(k)
- :
- 1.0);
+ const double J = data.JxW_values[k] / quadrature.weight(k);
data.shape_gradients[first+d][k] = sign_change[i] *
shape_grads2[k][d] / J;
}
{
switch (mapping_type)
{
- case independent:
- case independent_on_cartesian:
+ case mapping_none:
{
for (unsigned int k=0; k<n_q_points; ++k)
for (unsigned int d=0; d<dim; ++d)
break;
}
- case covariant:
- case contravariant:
+ case mapping_covariant:
+ case mapping_contravariant:
{
// Use auxiliary vector
// for transformation
std::vector<Tensor<1,dim> > shape_values (n_q_points);
- if (mapping_type == covariant)
+ if (mapping_type == mapping_covariant)
mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
shape_values, mapping_data, mapping_covariant);
else
// face in
// fill_fe_face_values
const double
- J = (mapping_type == contravariant ?
+ J = (mapping_type == mapping_contravariant ?
data.cell_JxW_values[k] / quadrature.weight(k)
:
1.0);
break;
}
+ case mapping_piola:
+ if (true)
+ {
+ std::vector<Tensor<1,dim> > shape_values (n_q_points);
+ mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
+ shape_values, mapping_data, mapping_piola);
+ for (unsigned int k=0; k<n_q_points; ++k)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_values(first+d,k)
+ = sign_change[i] * shape_values[k][d];
+ break;
+ }
default:
Assert(false, ExcNotImplemented());
switch (mapping_type)
{
- case independent:
- case independent_on_cartesian:
+ case mapping_none:
{
mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
shape_grads1, mapping_data, mapping_covariant);
break;
}
- case covariant:
+ case mapping_covariant:
{
mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
shape_grads1, mapping_data, mapping_covariant);
data.shape_gradients[first+d][k] = shape_grads1[k][d];
break;
}
-
- case contravariant:
+ case mapping_piola:
+ case mapping_contravariant:
{
mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
shape_grads1, mapping_data, mapping_covariant);
// recompute
// determinant in the
// same way as above
- const double
- J = data.cell_JxW_values[k] / quadrature.weight(k);
+ const double J = data.cell_JxW_values[k] / quadrature.weight(k);
data.shape_gradients[first+d][k]
= sign_change[i] * shape_grads2[k][d] / J;
}
{
switch (mapping_type)
{
- case independent:
- case independent_on_cartesian:
+ case mapping_none:
{
for (unsigned int k=0; k<n_q_points; ++k)
for (unsigned int d=0; d<dim; ++d)
break;
}
- case covariant:
- case contravariant:
+ case mapping_covariant:
+ case mapping_contravariant:
{
// Use auxiliary vector for
// transformation
std::vector<Tensor<1,dim> > shape_values (n_q_points);
- if (mapping_type == covariant)
+ if (mapping_type == mapping_covariant)
mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
shape_values, mapping_data, mapping_covariant);
else
// face in
// fill_fe_face_values
const double
- J = (mapping_type == contravariant ?
+ J = (mapping_type == mapping_contravariant ?
data.cell_JxW_values[k] / quadrature.weight(k)
:
1.0);
break;
}
-
+ case mapping_piola:
+ if (true)
+ {
+ std::vector<Tensor<1,dim> > shape_values (n_q_points);
+ mapping.transform(make_slice(fe_data.shape_values[i], offset, n_q_points),
+ shape_values, mapping_data, mapping_piola);
+ for (unsigned int k=0; k<n_q_points; ++k)
+ for (unsigned int d=0; d<dim; ++d)
+ data.shape_values(first+d,k)
+ = sign_change[i] * shape_values[k][d];
+ break;
+ }
default:
Assert(false, ExcNotImplemented());
}
switch (mapping_type)
{
- case independent:
- case independent_on_cartesian:
+ case mapping_none:
{
mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
shape_grads1, mapping_data, mapping_covariant);
break;
}
- case covariant:
+ case mapping_covariant:
{
mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
shape_grads1, mapping_data, mapping_covariant);
data.shape_gradients[first+d][k] = shape_grads1[k][d];
break;
}
-
- case contravariant:
+
+ case mapping_piola:
+ case mapping_contravariant:
{
mapping.transform(make_slice(fe_data.shape_grads[i], offset, n_q_points),
shape_grads1, mapping_data, mapping_covariant);
// recompute
// determinant in the
// same way as above
- const double
- J = data.cell_JxW_values[k] / quadrature.weight(k);
+ const double J = data.cell_JxW_values[k] / quadrature.weight(k);
data.shape_gradients[first+d][k]
= sign_change[i] * shape_grads2[k][d] / J;
}
UpdateFlags
FE_PolyTensor<POLY,dim,spacedim>::update_once (const UpdateFlags flags) const
{
- Assert (mapping_type != no_mapping, ExcNotInitialized());
- const bool values_once = (mapping_type == independent);
+ const bool values_once = (mapping_type == mapping_none);
UpdateFlags out = update_default;
UpdateFlags
FE_PolyTensor<POLY,dim,spacedim>::update_each (const UpdateFlags flags) const
{
- Assert (mapping_type != no_mapping, ExcNotInitialized());
- const bool values_once = (mapping_type == independent);
+ const bool values_once = (mapping_type == mapping_none);
UpdateFlags out = update_default;
// $Id$
// Version: $Name$
//
-// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Assert (dim >= 2, ExcImpossibleInDim(dim));
const unsigned int n_dofs = this->dofs_per_cell;
- this->mapping_type = this->contravariant;
+ this->mapping_type = mapping_piola;
// First, initialize the
// generalized support points and
// quadrature weights, since they
if (flags & update_values)
out |= update_values | update_covariant_transformation
- | update_contravariant_transformation
+ | update_piola
| update_cell_JxW_values
| update_JxW_values;
if (flags & update_gradients)
out |= update_gradients | update_covariant_transformation
- | update_contravariant_transformation
+ | update_piola
| update_cell_JxW_values
| update_JxW_values;
if (flags & update_hessians)
// $Id$
// Version: $Name$
//
-// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Assert (dim >= 2, ExcImpossibleInDim(dim));
const unsigned int n_dofs = this->dofs_per_cell;
- this->mapping_type = this->independent_on_cartesian;
+ this->mapping_type = mapping_none;
// These must be done first, since
// they change the evaluation of
// basis functions
data->update_once = update_once(update_flags);
data->update_each = update_each(update_flags);
data->update_flags = data->update_once | data->update_each;
-
+
return data;
}
data->update_once = update_once(update_flags);
data->update_each = update_each(update_flags);
data->update_flags = data->update_once | data->update_each;
-
+
return data;
}
data->update_once = update_once(update_flags);
data->update_each = update_each(update_flags);
data->update_flags = data->update_once | data->update_each;
-
+
return data;
}
// equal and are the product of the
// local lengths in each coordinate
// direction
- if (data.current_update_flags() & update_JxW_values)
+ if (data.current_update_flags() & (update_JxW_values | update_volume_elements))
{
double J = data.length[0];
for (unsigned int d=1;d<dim;++d)
J *= data.length[d];
- for (unsigned int i=0; i<JxW_values.size();++i)
- JxW_values[i] = J * q.weight(i);
+ data.volume_element = J;
+ if (data.current_update_flags() & update_JxW_values)
+ for (unsigned int i=0; i<JxW_values.size();++i)
+ JxW_values[i] = J * q.weight(i);
}
// "compute" Jacobian at the quadrature
// points, which are all the same
for (unsigned int i=0; i<boundary_forms.size();++i)
boundary_forms[i] = J * normal_vectors[i];
- //TODO: Implement cell_JxW_values
+ if (data.current_update_flags() & update_volume_elements)
+ {
+ J = data.length[0];
+ for (unsigned int d=1;d<dim;++d)
+ J *= data.length[d];
+ data.volume_element = J;
+ }
}
for (unsigned int i=0; i<boundary_forms.size();++i)
boundary_forms[i] = J * normal_vectors[i];
- //TODO: Implement cell_JxW_values
+ if (data.current_update_flags() & update_volume_elements)
+ {
+ J = data.length[0];
+ for (unsigned int d=1;d<dim;++d)
+ J *= data.length[d];
+ data.volume_element = J;
+ }
}
output[i][d] = input[i][d]*data.length[d];
return;
}
+ case mapping_piola:
+ if (true)
+ {
+ 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)
+ output[i][d] = input[i][d] * data.length[d] / data.volume_element;
+ return;
+ }
default:
Assert(false, ExcNotImplemented());
}
| update_cell_JxW_values
| update_boundary_forms
| update_normal_vectors
+ | update_volume_elements
| update_jacobians
| update_jacobian_grads
| update_inverse_jacobians));
if (flags & update_contravariant_transformation)
data.contravariant.resize(n_original_q_points);
+ if (flags & update_volume_elements)
+ data.volume_elements.resize(n_original_q_points);
+
if (flags & update_jacobian_grads)
data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
data.contravariant[point][i][j] += data_derv[j] * supp_pts[i];
}
+ if (update_flags & update_volume_elements)
+ for (unsigned int point=0; point<n_q_points; ++point)
+ data.volume_elements[point] = determinant(data.contravariant[point]);
// invert contravariant for
// covariant transformation
}
return;
}
+ case mapping_piola:
+ if (true)
+ {
+ 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)
+ auxiliary[d] = input[i][d] / data.volume_elements[i];
+ contract (output[i], data.contravariant[i], auxiliary);
+ }
+ return;
+ }
default:
Assert(false, ExcNotImplemented());
}