* transformation.
*/
update_jacobian_grads = 0x0100,
+ //! Volume element
+ /**
+ * Compute the inverse
+ * Jacobian of the
+ * transformation from the
+ * reference cell to the real
+ * cell.
+ */
+ update_inverse_jacobians = 0x0200,
//! Covariant transformation
/**
* Compute all values the
* mappings like
* MappingCartesian this may be
* simpler than
- * #update_jacobians.
+ * #update_inverse_jacobians.
*/
- update_covariant_transformation = 0x0200,
+ update_covariant_transformation = 0x0400,
//! Contravariant transformation
/**
* Compute all values the
* simpler than
* #update_jacobians.
*/
- update_contravariant_transformation = 0x0400,
+ update_contravariant_transformation = 0x0800,
//! Shape function values of transformation
/**
* Compute the shape function
* values of the transformation
* defined by the Mapping.
*/
- update_transformation_values = 0x0800,
+ update_transformation_values = 0x1000,
//! Shape function gradients of transformation
/**
* Compute the shape function
* transformation defined by
* the Mapping.
*/
- update_transformation_gradients = 0x1000,
+ update_transformation_gradients = 0x2000,
/**
* Compute the JxW values
* is used in conjunction with
* H_div subspaces like RT and ABF.
*/
- update_cell_JxW_values = 0x2000,
+ update_cell_JxW_values = 0x4000,
/**
* Update the location of the
* mapped generalized support
* points of the element.
*/
- update_support_points = 0x4000,
+ update_support_points = 0x8000,
/**
* Update the Jacobian of the
* mapping in generalized
* support points.
*/
- update_support_jacobians = 0x8000,
+ update_support_jacobians = 0x10000,
/**
* Update the inverse Jacobian
* of the mapping in
* generalized support points.
*/
- update_support_inverse_jacobians = 0x10000,
+ update_support_inverse_jacobians = 0x20000,
/**
* @deprecated Update
* quadrature points
* quadrature points.
*/
std::vector<Tensor<2,dim> > jacobians;
-
+
/**
* Array of the derivatives of the Jacobian
* matrices at the quadrature points.
*/
std::vector<Tensor<3,dim> > jacobian_grads;
-
+
+ /**
+ * Array of the inverse Jacobian matrices
+ * at the quadrature points.
+ */
+ std::vector<Tensor<2,dim> > inverse_jacobians;
+
/**
* Store an array of weights
* times the Jacobi determinant
* jacobian_grads().
*/
const std::vector<Tensor<3,dim> > & get_jacobian_grads () const;
-
+
+ /**
+ * Return the inverse Jacobian of the
+ * transformation at the specified
+ * quadrature point, i.e.
+ * $J_{ij}=d\hat x_i/dx_j$
+ */
+ const Tensor<2,dim> & inverse_jacobian (const unsigned int quadrature_point) const;
+
+ /**
+ * Pointer to the array holding
+ * the values returned by
+ * inverse_jacobian().
+ */
+ const std::vector<Tensor<2,dim> > & get_inverse_jacobians () const;
+
/**
* Constant reference to the
* selected mapping object.
+template <int dim>
+inline
+const std::vector<Tensor<2,dim> >&
+FEValuesBase<dim>::get_inverse_jacobians () const
+{
+ Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField());
+ return this->inverse_jacobians;
+}
+
+
+
template <int dim>
inline
const Point<dim> &
+template <int dim>
+inline
+const Tensor<2,dim> &
+FEValuesBase<dim>::inverse_jacobian (const unsigned int i) const
+{
+ Assert (this->update_flags & update_inverse_jacobians, ExcAccessToUninitializedField());
+ Assert (i<this->inverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size()));
+
+ return this->inverse_jacobians[i];
+}
+
+
+
template <int dim>
template <class InputVector>
inline
* matrices needed to transform
* vector-valued functions,
* namely
- * @p jacobians
- * and the derivatives
- * @p jacobian_grads.
+ * @p jacobians,
+ * the derivatives
+ * @p jacobian_grads,
+ * and the inverse operations in
+ * @p inverse_jacobians.
*/
virtual void
fill_fe_values (const typename Triangulation<dim>::cell_iterator &cell,
std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values,
std::vector<Tensor<2,dim> > &jacobians,
- std::vector<Tensor<3,dim> > &jacobian_grads) const = 0;
+ std::vector<Tensor<3,dim> > &jacobian_grads,
+ std::vector<Tensor<2,dim> > &inverse_jacobians) const = 0;
/**
* Performs the same as @p fill_fe_values
std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values,
std::vector<Tensor<2,dim> > &jacobians,
- std::vector<Tensor<3,dim> > &jacobian_grads) const ;
+ std::vector<Tensor<3,dim> > &jacobian_grads,
+ std::vector<Tensor<2,dim> > &inverse_jacobians) const ;
virtual void
fill_fe_face_values (const typename Triangulation<dim>::cell_iterator &cell,
typename std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values,
std::vector<Tensor<2,dim> > &jacobians,
- std::vector<Tensor<3,dim> > &jacobian_grads) const ;
+ std::vector<Tensor<3,dim> > &jacobian_grads,
+ std::vector<Tensor<2,dim> > &inverse_jacobians) const ;
/**
* Implementation of the interface in
typename std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values,
std::vector<Tensor<2,dim> > &jacobians,
- std::vector<Tensor<3,dim> > &jacobian_grads) const ;
+ std::vector<Tensor<3,dim> > &jacobian_grads,
+ std::vector<Tensor<2,dim> > &inverse_jacobians) const ;
/**
* Implementation of the interface in
if (flags & update_jacobian_grads)
this->jacobian_grads.resize(n_quadrature_points);
+ if (flags & update_inverse_jacobians)
+ this->inverse_jacobians.resize(n_quadrature_points);
+
if (flags & update_boundary_forms)
this->boundary_forms.resize(n_quadrature_points);
this->quadrature_points,
this->JxW_values,
this->jacobians,
- this->jacobian_grads);
+ this->jacobian_grads,
+ this->inverse_jacobians);
this->get_fe().fill_fe_values(this->get_mapping(),
*this->present_cell,
std::vector<Point<dim> >& quadrature_points,
std::vector<double>& JxW_values,
std::vector<Tensor<2,dim> >& jacobians,
- std::vector<Tensor<3,dim> >& jacobian_grads) const
+ std::vector<Tensor<3,dim> >& jacobian_grads,
+ std::vector<Tensor<2,dim> >& inverse_jacobians) const
{
// convert data object to internal
// data for this class. fails with
if (data.current_update_flags() & update_jacobian_grads)
for (unsigned int i=0; i<jacobian_grads.size();++i)
jacobian_grads[i]=Tensor<3,dim>();
+ // "compute" inverse Jacobian at the
+ // quadrature points, which are all
+ // the same
+ if (data.current_update_flags() & update_inverse_jacobians)
+ for (unsigned int i=0; i<inverse_jacobians.size();++i)
+ {
+ inverse_jacobians[i]=Tensor<2,dim>();
+ for (unsigned int j=0; j<dim; ++j)
+ inverse_jacobians[j][j]=1./data.length[j];
+ }
}
std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values,
std::vector<Tensor<2,dim> > &jacobians,
- std::vector<Tensor<3,dim> > &jacobian_grads) const
+ std::vector<Tensor<3,dim> > &jacobian_grads,
+ std::vector<Tensor<2,dim> > &inverse_jacobians) const
{
// convert data object to internal
// data for this class. fails with
MappingQ1<dim>::fill_fe_values(cell, q, *p_data,
quadrature_points, JxW_values,
- jacobians, jacobian_grads);
+ jacobians, jacobian_grads, inverse_jacobians);
}
| update_boundary_forms
| update_normal_vectors
| update_jacobians
- | update_jacobian_grads))
+ | update_jacobian_grads
+ | update_inverse_jacobians))
out |= update_transformation_gradients;
return out;
| update_boundary_forms
| update_normal_vectors
| update_jacobians
- | update_jacobian_grads));
+ | update_jacobian_grads
+ | update_inverse_jacobians));
// add a few flags. note that some
// flags appear in both conditions
// operations. this leads to some
// circular logic. the only way to
// treat this is to iterate. since
- // there are 3 if-clauses in the
+ // there are 4 if-clauses in the
// loop, it will take at most 3
// iterations to converge. do them:
- for (unsigned int i=0; i<3; ++i)
+ for (unsigned int i=0; i<4; ++i)
{
// The following is a little incorrect:
// If not applied on a face,
| update_boundary_forms
| update_normal_vectors))
out |= update_contravariant_transformation;
+
+ if (out & (update_inverse_jacobians))
+ out |= update_covariant_transformation;
// The contravariant transformation
// is a Piola transformation, which
template <int dim>
void
MappingQ1<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator &cell,
- const Quadrature<dim> &q,
- typename Mapping<dim>::InternalDataBase &mapping_data,
+ const Quadrature<dim> &q,
+ typename Mapping<dim>::InternalDataBase &mapping_data,
std::vector<Point<dim> > &quadrature_points,
std::vector<double> &JxW_values,
std::vector<Tensor<2,dim> > &jacobians,
- std::vector<Tensor<3,dim> > &jacobian_grads) const
+ std::vector<Tensor<3,dim> > &jacobian_grads,
+ std::vector<Tensor<2,dim> > &inverse_jacobians) const
{
InternalData *data_ptr = dynamic_cast<InternalData *> (&mapping_data);
Assert(data_ptr!=0, ExcInternalError());
+= (data.second_derivative(point+DataSetDescriptor::cell (), k)[j][l]
*
data.mapping_support_points[k][i]);
+ }
+ // 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));
+ for (unsigned int point=0; point<n_q_points; ++point)
+ inverse_jacobians[point]
+ = transpose(data.covariant[point]);
}
}