#include <cmath>
#include <base/point.h>
#include <base/subscriptor.h>
-#include <base/vector2d.h>
#include <grid/tria.h>
#include <dofs/dof_handler.h>
#include <fe/fe_update_flags.h>
};
/**
- * Tranform a field of covariant
- * vectors. The vector spanned by
- * @p{begin} and @p{end} must
+ * Transform a field of covariant
+ * vectors. The covariant
+ * transformation multiplies a
+ * vector from the right with the
+ * inverse of the Jacobian matrix
+ * of the transformation from
+ * unit to real space
+ * cell. Alternatively, this is
+ * equivalent to a multiplication
+ * from the left with the
+ * transpose if the inverse
+ * matrix.
+ *
+ * The range of vectors spanned
+ * by @p{begin} and @p{end} must
* have as many elements as there
* are quadrature points (not
* tested inside the function).
+ * Also note the different
+ * convention for parameters
+ * compared to the standard C++
+ * @p{transform} function: here,
+ * first destination, then source
+ * are specified, and the number
+ * of elements is determined by a
+ * range of destination
+ * vectors. This convention is
+ * due to the way the function is
+ * usually called.
*
* The vector @p{src} must
* contain at least as many
- * elements.
+ * elements as there are
+ * quadrature points.
*/
virtual
void
transform_covariant (Tensor<1,dim> *begin,
Tensor<1,dim> *end,
const Tensor<1,dim> *src,
+ const InternalDataBase &internal) const = 0;
+
+ /**
+ * Transform a set of matrices
+ * covariantly, i.e. multiply
+ * each function in the input
+ * range by the Jacobian matrices
+ * at the different quadrature
+ * points from the left.
+ *
+ * The same applies as to the
+ * function above regarding input
+ * and output ranges.
+ */
+ virtual
+ void
+ transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
const InternalDataBase &internal) const = 0;
/**
- * Tranform a field of
+ * Transform a field of
* contravariant vectors. The
- * vector spanned by @p{begin}
- * and @p{end} must have as many
- * elements as there are
- * quadrature points (not tested
- * inside the function).
+ * contravariant transformation
+ * multiplies a vector from the
+ * left with the Jacobian matrix
+ * of the transformation from
+ * unit to real space cell.
+ *
+ * The range of vectors spanned
+ * by @p{begin} and @p{end} must
+ * have as many elements as there
+ * are quadrature points (not
+ * tested inside the function).
+ * Also note the different
+ * convention for parameters
+ * compared to the standard C++
+ * @p{transform} function: here,
+ * first destination, then source
+ * are specified, and the number
+ * of elements is determined by a
+ * range of destination
+ * vectors. This convention is
+ * due to the way the function is
+ * usually called.
*
* The vector @p{src} must
* contain at least as many
- * elements.
+ * elements as there are
+ * quadrature points.
*/
virtual
void
const Tensor<1,dim> *src,
const InternalDataBase &internal) const = 0;
+ /**
+ * Transform a set of matrices
+ * contravariantly, i.e. multiply
+ * each function in the input
+ * range by the inverse Jacobian
+ * matrices at the different
+ * quadrature points from the
+ * right. Note that here it is no
+ * more equivalent to
+ * multiplication with the
+ * transpose of the inverse
+ * matrices from the left, unless
+ * the matrices to be multiplied
+ * with are symmetric.
+ *
+ * The same applies as to the
+ * function above regarding input
+ * and output ranges.
+ */
+ virtual
+ void
+ transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const InternalDataBase &internal) const = 0;
+
/**
* Indicate fields to be updated
* in the constructor of
fill_fe_values (const typename DoFHandler<dim>::cell_iterator &cell,
const Quadrature<dim> &quadrature,
InternalDataBase &internal,
- std::vector<Point<dim> > &quadrature_points,
+ std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values) const = 0;
/**
transform_covariant (Tensor<1,dim> *begin,
Tensor<1,dim> *end,
const Tensor<1,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &internal) const;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
const typename Mapping<dim>::InternalDataBase &internal) const;
/**
transform_contravariant (Tensor<1,dim> *begin,
Tensor<1,dim> *end,
const Tensor<1,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &internal) const;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
const typename Mapping<dim>::InternalDataBase &internal) const;
/**
*
* Filled once.
*/
- vector2d<Tensor<1,dim> > unit_tangentials;
+ Table<2,Tensor<1,dim> > unit_tangentials;
/**
* Auxiliary vectors for internal use.
*/
- vector2d<Tensor<1,dim> > aux;
+ Table<2,Tensor<1,dim> > aux;
};
/**
#include <base/config.h>
-#include <base/vector2d.h>
+#include <base/table.h>
#include <fe/mapping_q1.h>
template <int dim> class TensorProductPolynomials;
const Point<dim> &p) const;
/**
- * Implementation of the interface in
- * @ref{Mapping}.
+ * Implementation of the
+ * interface in @ref{Mapping}.
*/
virtual void
transform_covariant (Tensor<1,dim> *begin,
const typename Mapping<dim>::InternalDataBase &internal) const;
/**
- * Implementation of the interface in
- * @ref{Mapping}.
+ * Implementation of the
+ * interface in @ref{Mapping}.
+ */
+ virtual void
+ transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &internal) const;
+
+ /**
+ * Implementation of the
+ * interface in @ref{Mapping}.
*/
virtual void
transform_contravariant (Tensor<1,dim> *begin,
const Tensor<1,dim> *src,
const typename Mapping<dim>::InternalDataBase &internal) const;
+ /**
+ * Implementation of the
+ * interface in @ref{Mapping}.
+ */
+ virtual void
+ transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &internal) const;
+
/**
* Return the degree of the
* mapping, i.e. the value which
* this vector is computed.
*/
void
- set_laplace_on_quad_vector(vector2d<double> &loqvs) const;
+ set_laplace_on_quad_vector(Table<2,double> &loqvs) const;
/**
* This function is needed by the
* @p{degree>2} this vector is
* computed.
*/
- void set_laplace_on_hex_vector(vector2d<double> &lohvs) const;
+ void set_laplace_on_hex_vector(Table<2,double> &lohvs) const;
/**
* Computes the @p{laplace_on_quad(hex)_vector}.
* functions if the data is not
* yet hardcoded.
*/
- void compute_laplace_vector(vector2d<double> &lvs) const;
+ void compute_laplace_vector(Table<2,double> &lvs) const;
/**
* Takes a
* @p{n_inner} computed inner
* points are appended.
*/
- void apply_laplace_vector(const vector2d<double> &lvs,
+ void apply_laplace_vector(const Table<2,double> &lvs,
std::vector<Point<dim> > &a) const;
/**
* unit_support_points on the
* boundary of the quad
*/
- vector2d<double> laplace_on_quad_vector;
+ Table<2,double> laplace_on_quad_vector;
/**
* Needed by the
* (for @p{dim==3}). Filled by the
* constructor.
*/
- vector2d<double> laplace_on_hex_vector;
+ Table<2,double> laplace_on_hex_vector;
/**
* Exception.
const std::vector<Point<1> > &unit_points,
MappingQ1<1>::InternalData &data) const;
template <> void MappingQ<1>::set_laplace_on_quad_vector(
- vector2d<double> &) const;
+ Table<2,double> &) const;
template <> void MappingQ<3>::set_laplace_on_hex_vector(
- vector2d<double> &lohvs) const;
+ Table<2,double> &lohvs) const;
template <> void MappingQ<1>::compute_laplace_vector(
- vector2d<double> &) const;
+ Table<2,double> &) const;
template <> void MappingQ<1>::add_line_support_points (
const Triangulation<1>::cell_iterator &,
std::vector<Point<1> > &) const;
#include <base/config.h>
+#include <base/table.h>
#include <cmath>
#include <fe/mapping.h>
* @ref{Mapping}.
*/
virtual void
+ transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &internal) const;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
transform_contravariant (Tensor<1,dim> *begin,
Tensor<1,dim> *end,
const Tensor<1,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &internal) const;
+
+ /**
+ * Implementation of the interface in
+ * @ref{Mapping}.
+ */
+ virtual void
+ transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
const typename Mapping<dim>::InternalDataBase &internal) const;
/**
/**
* Tensors of covariant
- * transformation.
+ * transformation at each of
+ * the quadrature points. The
+ * matrix stored is the
+ * inverse of the Jacobian
+ * matrix, which itself is
+ * stored in the
+ * @p{contravariant} field of
+ * this structure.
*
* Computed on each cell.
*/
/**
* Tensors of covariant
- * transformation.
+ * transformation at each of
+ * the quadrature points. The
+ * contravariant matrix is
+ * the Jacobian of the
+ * transformation,
+ * i.e. $J_ij=dx_i/d\hat x_j$.
*
* Computed on each cell.
*/
*
* Filled once.
*/
- vector2d<Tensor<1,dim> > unit_tangentials;
+ Table<2,Tensor<1,dim> > unit_tangentials;
/**
* Auxuliary vectors for internal use.
*/
- vector2d<Tensor<1,dim> > aux;
+ Table<2,Tensor<1,dim> > aux;
/**
* Stores the support points of
* can be initialized by
* @p{DoFObjectAccessor::set_dof_values()}.
*/
- MappingQ1Eulerian ( const Vector<double> &euler_transform_vectors,
- const DoFHandler<dim> &shiftmap_dof_handler);
+ MappingQ1Eulerian (const Vector<double> &euler_transform_vectors,
+ const DoFHandler<dim> &shiftmap_dof_handler);
/**
}
+
+template <int dim>
+void
+MappingCartesian<dim>::transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+ const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
+
+ Assert (data.update_flags & update_covariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+ // simply scale by inverse Jacobian
+ // (which is diagonal here)
+ while (begin!=end)
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int p=0; p<dim; ++p)
+ (*begin)[d][p] = (*src)[d][p] / data.length[p];
+ begin++;
+ src++;
+ }
+}
+
+
+
template <int dim>
void
MappingCartesian<dim>::transform_contravariant (Tensor<1,dim> *begin,
}
+
+template <int dim>
+void
+MappingCartesian<dim>::transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+ // convert data object to internal
+ // data for this class. fails with
+ // an exception if that is not
+ // possible
+ const InternalData &data = dynamic_cast<const InternalData&> (mapping_data);
+
+ Assert (data.update_flags & update_contravariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+ // simply scale by Jacobian
+ // (which is diagonal here)
+ while (begin!=end)
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int p=0; p<dim; ++p)
+ (*begin)[d][p] = data.length[d] * (*src)[d][p];
+ begin++;
+ src++;
+ }
+}
+
+
+
template <int dim>
Point<dim> MappingCartesian<dim>::transform_unit_to_real_cell (
const typename Triangulation<dim>::cell_iterator cell,
// cells are scaled linearly
template<>
MappingQ<1>::MappingQ (const unsigned int):
- laplace_on_quad_vector(0),
- laplace_on_hex_vector(0),
degree(1),
n_inner(0),
n_outer(0),
template<int dim>
-MappingQ<dim>::MappingQ (const unsigned int p):
- laplace_on_quad_vector(0),
- laplace_on_hex_vector(0),
+MappingQ<dim>::MappingQ (const unsigned int p)
+ :
degree(p),
n_inner(power(degree-1, dim)),
n_outer((dim==2) ? 4+4*(degree-1)
template <>
void
-MappingQ<1>::set_laplace_on_quad_vector(vector2d<double> &) const
+MappingQ<1>::set_laplace_on_quad_vector(Table<2,double> &) const
{
Assert(false, ExcInternalError());
}
template <int dim>
void
-MappingQ<dim>::set_laplace_on_quad_vector(vector2d<double> &loqvs) const
+MappingQ<dim>::set_laplace_on_quad_vector(Table<2,double> &loqvs) const
{
Assert(degree>1, ExcInternalError());
const unsigned int n_inner_2d=(degree-1)*(degree-1);
template <>
void
-MappingQ<3>::set_laplace_on_hex_vector(vector2d<double> &lohvs) const
+MappingQ<3>::set_laplace_on_hex_vector(Table<2,double> &lohvs) const
{
Assert(degree>1, ExcInternalError());
template <int dim>
void
-MappingQ<dim>::set_laplace_on_hex_vector(vector2d<double> &) const
+MappingQ<dim>::set_laplace_on_hex_vector(Table<2,double> &) const
{
Assert(false, ExcInternalError());
}
template <>
void
-MappingQ<1>::compute_laplace_vector(vector2d<double> &) const
+MappingQ<1>::compute_laplace_vector(Table<2,double> &) const
{
Assert(false, ExcInternalError());
}
template <int dim>
void
-MappingQ<dim>::compute_laplace_vector(vector2d<double> &lvs) const
+MappingQ<dim>::compute_laplace_vector(Table<2,double> &lvs) const
{
Assert(lvs.n_rows()==0, ExcInternalError());
Assert(dim==2 || dim==3, ExcNotImplemented());
template <int dim>
void
-MappingQ<dim>::apply_laplace_vector(const vector2d<double> &lvs,
+MappingQ<dim>::apply_laplace_vector(const Table<2,double> &lvs,
std::vector<Point<dim> > &a) const
{
Assert(lvs.n_rows()!=0, ExcLaplaceVectorNotSet(degree));
}
+
+template <int dim>
+void
+MappingQ<dim>::transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+ const typename MappingQ1<dim>::InternalData *q1_data =
+ dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
+ Assert(q1_data!=0, ExcInternalError());
+
+ typename std::vector<Tensor<2,dim> >::const_iterator tensor;
+
+ if (q1_data->is_mapping_q1_data)
+ tensor = q1_data->covariant.begin();
+ else
+ {
+ const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
+ Assert(data!=0, ExcInternalError());
+
+ if (data->use_mapping_q1_on_current_cell)
+ tensor = data->mapping_q1_data.covariant.begin();
+ else
+ tensor = data->covariant.begin();
+ }
+ while (begin!=end)
+ {
+ contract (*(begin++), *(src++), *(tensor++));
+ }
+}
+
+
+
template <int dim>
void
MappingQ<dim>::transform_contravariant (Tensor<1,dim> *begin,
}
+
+template <int dim>
+void
+MappingQ<dim>::transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+ const typename MappingQ1<dim>::InternalData *q1_data =
+ dynamic_cast<const typename MappingQ1<dim>::InternalData *> (&mapping_data);
+ Assert(q1_data!=0, ExcInternalError());
+
+ typename std::vector<Tensor<2,dim> >::const_iterator tensor;
+
+ if (q1_data->is_mapping_q1_data)
+ tensor = q1_data->contravariant.begin();
+ else
+ {
+ const InternalData *data = dynamic_cast<const InternalData *> (q1_data);
+ Assert(data!=0, ExcInternalError());
+
+ if (data->use_mapping_q1_on_current_cell)
+ tensor = data->mapping_q1_data.contravariant.begin();
+ else
+ tensor = data->contravariant.begin();
+ }
+ while (begin!=end)
+ {
+ contract (*(begin++), *(tensor++), *(src++));
+ }
+}
+
+
+
template <int dim>
Point<dim> MappingQ<dim>::transform_unit_to_real_cell (
const typename Triangulation<dim>::cell_iterator cell,
// ExcDimensionMismatch(covariant_grads.size(), npts));
}
- std::vector<Point<dim> > &a=data.mapping_support_points;
-
- // store all Lagrangian
- // support points in a
- if (a.size()==0
- || (&cell->get_triangulation() !=
- &data.cell_of_current_support_points->get_triangulation())
- || (cell!=data.cell_of_current_support_points))
+ // if necessary, recompute the
+ // support points of the
+ // transformation of this cell
+ if ((data.mapping_support_points.size() == 0)
+ ||
+ (&cell->get_triangulation() !=
+ &data.cell_of_current_support_points->get_triangulation())
+ ||
+ (cell != data.cell_of_current_support_points))
{
- compute_mapping_support_points(cell, a);
- data.cell_of_current_support_points=cell;
+ compute_mapping_support_points(cell, data.mapping_support_points);
+ data.cell_of_current_support_points = cell;
}
+
+ // first compute quadrature points
+ if (update_flags & update_q_points)
+ for (unsigned int point=0; point<npts; ++point)
+ for (unsigned int k=0; k<data.n_shape_functions; ++k)
+ quadrature_points[point]
+ += data.shape(point+offset,k) * data.mapping_support_points[k];
- for (unsigned int point=0; point<npts; ++point)
- {
- // First, compute function
- // values and derivatives
+ // then Jacobians
+ if (update_flags & update_contravariant_transformation)
+ for (unsigned int point=0; point<npts; ++point)
for (unsigned int k=0; k<data.n_shape_functions; ++k)
- {
- if (update_flags & update_q_points)
- {
- quadrature_points[point]
- += data.shape(point+offset,k) * a[k];
- }
- if (update_flags & update_contravariant_transformation)
- {
- for (unsigned int i=0; i<dim; ++i)
- for (unsigned int j=0; j<dim; ++j)
- data.contravariant[point][i][j]
- += data.derivative(point+offset, k)[j]
- * a[k][i];
- }
- }
-
- // Invert contravariant for
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ data.contravariant[point][i][j]
+ += (data.derivative(point+offset, k)[j]
+ *
+ data.mapping_support_points[k][i]);
+
+ // invert contravariant for
// covariant transformation
// matrices
- if (update_flags & update_covariant_transformation)
- data.covariant[point]
- = invert(data.contravariant[point]);
- }
+ if (update_flags & update_covariant_transformation)
+ for (unsigned int point=0; point<npts; ++point)
+ data.covariant[point] = invert(data.contravariant[point]);
+
data.first_cell = false;
}
// Assert (dst.size() == data.contravariant.size(),
// ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size()));
- typename std::vector<Tensor<2,dim> >::const_iterator tensor = data.covariant.begin();
+ typename std::vector<Tensor<2,dim> >::const_iterator
+ tensor = data.covariant.begin();
+
+ while (begin!=end)
+ {
+ contract (*(begin++), *(src++), *(tensor++));
+ }
+}
+
+
+
+template <int dim>
+void
+MappingQ1<dim>::transform_covariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+ const InternalData *data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
+ Assert(data_ptr!=0, ExcInternalError());
+ const InternalData &data=*data_ptr;
+
+ Assert (data.update_flags & update_covariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+//TODO: [GK] Can we do a similar assertion?
+// Assert (dst.size() == data.contravariant.size(),
+// ExcDimensionMismatch(dst.size() + src_offset, data.contravariant.size()));
+
+ typename std::vector<Tensor<2,dim> >::const_iterator
+ tensor = data.covariant.begin();
while (begin!=end)
{
}
+
template <int dim>
void
MappingQ1<dim>::transform_contravariant (Tensor<1,dim> *begin,
Assert (data.update_flags & update_contravariant_transformation,
typename FEValuesBase<dim>::ExcAccessToUninitializedField());
- typename std::vector<Tensor<2,dim> >::const_iterator tensor = data.contravariant.begin();
+ typename std::vector<Tensor<2,dim> >::const_iterator
+ tensor = data.contravariant.begin();
while (begin!=end)
{
}
}
+
+template <int dim>
+void
+MappingQ1<dim>::transform_contravariant (Tensor<2,dim> *begin,
+ Tensor<2,dim> *end,
+ const Tensor<2,dim> *src,
+ const typename Mapping<dim>::InternalDataBase &mapping_data) const
+{
+ const InternalData* data_ptr = dynamic_cast<const InternalData *> (&mapping_data);
+ Assert(data_ptr!=0, ExcInternalError());
+ const InternalData &data=*data_ptr;
+
+ Assert (data.update_flags & update_contravariant_transformation,
+ typename FEValuesBase<dim>::ExcAccessToUninitializedField());
+
+ typename std::vector<Tensor<2,dim> >::const_iterator
+ tensor = data.contravariant.begin();
+
+ while (begin!=end)
+ {
+ contract (*(begin++), *(tensor++), *(src++));
+ }
+}
+
+
+
template <int dim>
Point<dim> MappingQ1<dim>::transform_unit_to_real_cell (
const typename Triangulation<dim>::cell_iterator cell,