* non-zero components.
*/
std::vector<unsigned int> shape_function_to_row_table;
-
+
+ /**
+ * Vector containing the permutation of
+ * shape functions necessary if the faces
+ * of a cell have the wrong
+ * face_orientation. This is computed only
+ * once. Actually, this does not contain
+ * the permutation itself but rather the
+ * shift of indices needed to calculate the
+ * permutation.
+ */
+ std::vector<int> shift_in_face_shape_functions;
+
+ /**
+ * Vector containing the permutation of
+ * shape functions due to faces with
+ * non-standard face_orientation on a given
+ * cell (recomputed on each cell.) If all
+ * faces are oriented according to the
+ * standard, this is the identity mapping.
+ */
+ std::vector<unsigned int> permutated_shape_functions;
+
+ /**
+ * Bool flag indicating the need to update
+ * the @p permutated_shape_functions vector
+ * on each cell. This is only necessary in
+ * 3d and if the finite element has
+ * shape_functions on the face.
+ */
+ bool update_shape_function_permutation;
+
/**
* Original update flags handed
* to the constructor of
const unsigned int point_no,
const unsigned int component) const;
-
+ /**
+ * If shape functions belong to a face in
+ * 3D, they have to be permutated, if the
+ * face has non-standard face
+ * orientation. This functuion takes an
+ * index of a shape function (on a standard
+ * cell) and returns the corresponding
+ * shape function on the real cell.
+ */
+ unsigned int
+ shift_shape_function_index (const unsigned int i) const;
+
+
//@}
/// @name FunctionAccess Access to values of global finite element functions
//@{
*/
UpdateFlags compute_update_flags (const UpdateFlags update_flags) const;
+ /**
+ * Reinit the permutation of (face) shape
+ * functions to match the present cell.
+ */
+ void reinit();
+
private:
/**
* Copy constructor. Since
FEValuesBase<dim>::shape_value (const unsigned int i,
const unsigned int j) const
{
+ const unsigned int I=shift_shape_function_index(i);
+
Assert (this->update_flags & update_values,
ExcAccessToUninitializedField());
- Assert (fe->is_primitive (i),
- ExcShapeFunctionNotPrimitive(i));
+ Assert (fe->is_primitive (I),
+ ExcShapeFunctionNotPrimitive(I));
// if the entire FE is primitive,
// then we can take a short-cut:
if (fe->is_primitive())
- return this->shape_values(i,j);
+ return this->shape_values(I,j);
else
// otherwise, use the mapping
// between shape function numbers
// question to which vector
// component the call of this
// function refers
- return this->shape_values(this->shape_function_to_row_table[i], j);
+ return this->shape_values(this->shape_function_to_row_table[I], j);
}
const unsigned int j,
const unsigned int component) const
{
+ const unsigned int I=shift_shape_function_index(i);
+
Assert (this->update_flags & update_values,
ExcAccessToUninitializedField());
Assert (component < fe->n_components(),
// system_to_component_table only
// works if the shape function is
// primitive):
- if (fe->is_primitive(i))
+ if (fe->is_primitive(I))
{
- if (component == fe->system_to_component_index(i).first)
- return this->shape_values(this->shape_function_to_row_table[i],j);
+ if (component == fe->system_to_component_index(I).first)
+ return this->shape_values(this->shape_function_to_row_table[I],j);
else
return 0;
}
// whether the shape function
// is non-zero at all within
// this component:
- if (fe->get_nonzero_components(i)[component] == false)
+ if (fe->get_nonzero_components(I)[component] == false)
return 0.;
// count how many non-zero
// shape function in the arrays
// we index presently:
const unsigned int
- row = (this->shape_function_to_row_table[i]
+ row = (this->shape_function_to_row_table[I]
+
- std::count (fe->get_nonzero_components(i).begin(),
- fe->get_nonzero_components(i).begin()+component,
+ std::count (fe->get_nonzero_components(I).begin(),
+ fe->get_nonzero_components(I).begin()+component,
true));
return this->shape_values(row, j);
};
FEValuesBase<dim>::shape_grad (const unsigned int i,
const unsigned int j) const
{
+ const unsigned int I=shift_shape_function_index(i);
+
Assert (this->update_flags & update_gradients,
ExcAccessToUninitializedField());
- Assert (fe->is_primitive (i),
- ExcShapeFunctionNotPrimitive(i));
+ Assert (fe->is_primitive (I),
+ ExcShapeFunctionNotPrimitive(I));
Assert (i<this->shape_gradients.size(),
- ExcIndexRange (i, 0, this->shape_gradients.size()));
+ ExcIndexRange (I, 0, this->shape_gradients.size()));
Assert (j<this->shape_gradients[0].size(),
ExcIndexRange (j, 0, this->shape_gradients[0].size()));
// if the entire FE is primitive,
// then we can take a short-cut:
if (fe->is_primitive())
- return this->shape_gradients[i][j];
+ return this->shape_gradients[I][j];
else
// otherwise, use the mapping
// between shape function numbers
// question to which vector
// component the call of this
// function refers
- return this->shape_gradients[this->shape_function_to_row_table[i]][j];
+ return this->shape_gradients[this->shape_function_to_row_table[I]][j];
}
const unsigned int j,
const unsigned int component) const
{
+ const unsigned int I=shift_shape_function_index(i);
+
Assert (this->update_flags & update_gradients,
ExcAccessToUninitializedField());
Assert (component < fe->n_components(),
// system_to_component_table only
// works if the shape function is
// primitive):
- if (fe->is_primitive(i))
+ if (fe->is_primitive(I))
{
- if (component == fe->system_to_component_index(i).first)
- return this->shape_gradients[this->shape_function_to_row_table[i]][j];
+ if (component == fe->system_to_component_index(I).first)
+ return this->shape_gradients[this->shape_function_to_row_table[I]][j];
else
return Tensor<1,dim>();
}
// whether the shape function
// is non-zero at all within
// this component:
- if (fe->get_nonzero_components(i)[component] == false)
+ if (fe->get_nonzero_components(I)[component] == false)
return Tensor<1,dim>();
// count how many non-zero
// shape function in the arrays
// we index presently:
const unsigned int
- row = (this->shape_function_to_row_table[i]
+ row = (this->shape_function_to_row_table[I]
+
- std::count (fe->get_nonzero_components(i).begin(),
- fe->get_nonzero_components(i).begin()+component,
+ std::count (fe->get_nonzero_components(I).begin(),
+ fe->get_nonzero_components(I).begin()+component,
true));
return this->shape_gradients[row][j];
};
FEValuesBase<dim>::shape_2nd_derivative (const unsigned int i,
const unsigned int j) const
{
+ const unsigned int I=shift_shape_function_index(i);
+
Assert (this->update_flags & update_second_derivatives,
ExcAccessToUninitializedField());
- Assert (fe->is_primitive (i),
+ Assert (fe->is_primitive (I),
ExcShapeFunctionNotPrimitive(i));
- Assert (i<this->shape_2nd_derivatives.size(),
- ExcIndexRange (i, 0, this->shape_2nd_derivatives.size()));
+ Assert (I<this->shape_2nd_derivatives.size(),
+ ExcIndexRange (I, 0, this->shape_2nd_derivatives.size()));
Assert (j<this->shape_2nd_derivatives[0].size(),
ExcIndexRange (j, 0, this->shape_2nd_derivatives[0].size()));
// if the entire FE is primitive,
// then we can take a short-cut:
if (fe->is_primitive())
- return this->shape_2nd_derivatives[i][j];
+ return this->shape_2nd_derivatives[I][j];
else
// otherwise, use the mapping
// between shape function numbers
// question to which vector
// component the call of this
// function refers
- return this->shape_2nd_derivatives[this->shape_function_to_row_table[i]][j];
+ return this->shape_2nd_derivatives[this->shape_function_to_row_table[I]][j];
}
const unsigned int j,
const unsigned int component) const
{
+ const unsigned int I=shift_shape_function_index(i);
+
Assert (this->update_flags & update_second_derivatives,
ExcAccessToUninitializedField());
Assert (component < fe->n_components(),
// system_to_component_table only
// works if the shape function is
// primitive):
- if (fe->is_primitive(i))
+ if (fe->is_primitive(I))
{
- if (component == fe->system_to_component_index(i).first)
- return this->shape_2nd_derivatives[this->shape_function_to_row_table[i]][j];
+ if (component == fe->system_to_component_index(I).first)
+ return this->shape_2nd_derivatives[this->shape_function_to_row_table[I]][j];
else
return Tensor<2,dim>();
}
// shape function in the arrays
// we index presently:
const unsigned int
- row = (this->shape_function_to_row_table[i]
+ row = (this->shape_function_to_row_table[I]
+
- std::count (fe->get_nonzero_components(i).begin(),
- fe->get_nonzero_components(i).begin()+component,
+ std::count (fe->get_nonzero_components(I).begin(),
+ fe->get_nonzero_components(I).begin()+component,
true));
return this->shape_2nd_derivatives[row][j];
};
}
+
+template <int dim>
+inline
+unsigned int
+FEValuesBase<dim>::shift_shape_function_index (const unsigned int i) const
+{
+ // standard implementation for 1D and 2D
+ Assert(i<fe->dofs_per_cell, ExcInternalError());
+ return i;
+}
+
+template <>
+inline
+unsigned int
+FEValuesBase<3>::shift_shape_function_index (const unsigned int i) const
+{
+ Assert(i<fe->dofs_per_cell, ExcInternalError());
+ return this->permutated_shape_functions[i];
+}
+
+
+
/*------------------------ Inline functions: FEValues ----------------------------*/
return this->boundary_forms[i];
}
+
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
if (flags & update_cell_JxW_values)
this->cell_JxW_values.resize(n_quadrature_points);
+
+ // initialize the permutation fields, if they
+ // are needed
+ if (dim==3)
+ {
+ permutated_shape_functions.resize(fe.dofs_per_cell);
+ // initialize cell permutation mapping
+ // with identity
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ this->permutated_shape_functions[i]=i;
+
+ if (fe.dofs_per_quad>0 &&
+ (flags & update_values ||
+ flags & update_gradients ||
+ flags & update_second_derivatives))
+ {
+ // ask the fe to fill the vector of
+ // shifts
+ fe.get_face_shape_function_shifts(shift_in_face_shape_functions);
+ Assert (shift_in_face_shape_functions.size()==fe.dofs_per_quad,
+ ExcInternalError());
+ update_shape_function_permutation=true;
+ }
+ else
+ update_shape_function_permutation=false;
+ }
+ else
+ update_shape_function_permutation=false;
}
}
+#if deal_II_dimension <3
+template <int dim>
+void FEValuesBase<dim>::reinit ()
+{
+ // do nothing in 1D and 2D
+}
+#else
+
+template <>
+void FEValuesBase<3>::reinit ()
+{
+ // in 3D reinit the permutation of face shape
+ // functions, if necessary
+ if (update_shape_function_permutation)
+ {
+ Assert (this->shift_in_face_shape_functions.size()==this->fe->dofs_per_quad,
+ ExcInternalError());
+ Triangulation<3>::cell_iterator thiscell = *(this->present_cell);
+ unsigned int offset=fe->first_quad_index;
+ for (unsigned int face_no=0; face_no<GeometryInfo<3>::faces_per_cell; ++face_no)
+ {
+ if (thiscell->face_orientation(face_no))
+ for (unsigned int i=0; i<fe->dofs_per_quad; ++i)
+ this->permutated_shape_functions[offset+i]=offset+i;
+ else
+ for (unsigned int i=0; i<fe->dofs_per_quad; ++i)
+ this->permutated_shape_functions[offset+i]=offset+i+this->shift_in_face_shape_functions[i];
+
+ offset+=this->fe->dofs_per_quad;
+ }
+ Assert (offset-fe->first_quad_index==GeometryInfo<3>::faces_per_cell*this->fe->dofs_per_quad,
+ ExcInternalError());
+ }
+}
+#endif
+
+
template <int dim>
template <class InputVector, typename number>
return (MemoryConsumption::memory_consumption (this->shape_values) +
MemoryConsumption::memory_consumption (this->shape_gradients) +
MemoryConsumption::memory_consumption (this->shape_2nd_derivatives) +
+ MemoryConsumption::memory_consumption (this->permutated_shape_functions) +
+ MemoryConsumption::memory_consumption (this->shift_in_face_shape_functions) +
MemoryConsumption::memory_consumption (this->JxW_values) +
MemoryConsumption::memory_consumption (this->quadrature_points) +
MemoryConsumption::memory_consumption (this->normal_vectors) +
MemoryConsumption::memory_consumption (this->boundary_forms) +
MemoryConsumption::memory_consumption (this->cell_JxW_values) +
sizeof(this->update_flags) +
+ sizeof(this->update_shape_function_permutation) +
MemoryConsumption::memory_consumption (n_quadrature_points) +
MemoryConsumption::memory_consumption (dofs_per_cell) +
MemoryConsumption::memory_consumption (mapping) +
this->fe_data->clear_first_cell ();
this->mapping_data->clear_first_cell ();
+
+ FEValuesBase<dim>::reinit();
}
this->fe_data->clear_first_cell ();
this->mapping_data->clear_first_cell ();
+ FEValuesBase<dim>::reinit();
}
this->fe_data->clear_first_cell ();
this->mapping_data->clear_first_cell ();
+ FEValuesBase<dim>::reinit();
}