#include <deal.II/base/aligned_vector.h>
#include <deal.II/fe/fe.h>
-#include <deal.II/matrix_free/helper_functions.h>
-
DEAL_II_NAMESPACE_OPEN
AlignedVector<VectorizedArray<Number> > shape_hessians_collocation_eo;
/**
- * Stores the indices from cell DoFs to face DoFs. The rows go through
- * the <tt>2*dim</tt> faces, and the columns the DoFs on the faces.
- */
- dealii::Table<2,unsigned int> face_indices;
-
- /**
- * Stores one-dimensional values of shape functions evaluated in zero
- * and one, i.e., on the one-dimensional faces. Not vectorized.
+ * Collects all data of 1D shape values evaluated at the point 0 and 1
+ * (the vertices) in one data structure. Sorting is first the values,
+ * then gradients, then second derivatives.
*/
- std::vector<Number> face_value[2];
-
- /**
- * Stores one-dimensional gradients of shape functions evaluated in zero
- * and one, i.e., on the one-dimensional faces. Not vectorized.
- */
- std::vector<Number> face_gradient[2];
+ AlignedVector<VectorizedArray<Number> > shape_data_on_face[2];
/**
* Stores one-dimensional values of shape functions on subface. Since
- * there are two subfaces, store two variants. Not vectorized.
+ * there are two subfaces, store two variants.
*/
- std::vector<Number> subface_value[2];
+ AlignedVector<VectorizedArray<Number> > values_within_subface[2];
/**
- * Non-vectorized version of shape values. Needed when evaluating face
- * info.
+ * Stores one-dimensional gradients of shape functions on subface. Since
+ * there are two subfaces, store two variants.
*/
- std::vector<Number> shape_values_number;
+ AlignedVector<VectorizedArray<Number> > gradients_within_subface[2];
/**
- * Non-vectorized version of shape gradients. Needed when evaluating
- * face info.
+ * Stores one-dimensional gradients of shape functions on subface. Since
+ * there are two subfaces, store two variants.
*/
- std::vector<Number> shape_gradient_number;
+ AlignedVector<VectorizedArray<Number> > hessians_within_subface[2];
/**
* Renumbering from deal.II's numbering of cell degrees of freedom to
*/
unsigned int dofs_per_face;
+ /**
+ * Indicates whether the basis functions are nodal in 0 and 1, i.e., the
+ * end points of the unit cell.
+ */
+ bool nodal_at_cell_boundaries;
+
+ /**
+ * For nodal cells, we might get around by simply loading the indices to
+ * the degrees of freedom that act on a particular face, rather than the
+ * whole set of indices that is then interpolated down to the
+ * element. This array stores this indirect addressing.
+ *
+ * The first table index runs through the faces of a cell, and the
+ * second runs through the nodal degrees of freedom of the face, using
+ * @p dofs_per_face entries.
+ *
+ * @note This object is only filled in case @p nodal_at_cell_boundaries
+ * evaluates to @p true.
+ */
+ dealii::Table<2,unsigned int> face_to_cell_index_nodal;
+
+ /**
+ * For Hermite-type basis functions, the @p face_to_cell_index_nodal for
+ * the values on a face of the cell is used together with a respective
+ * slot in the derivative. In the lexicographic ordering, this index is
+ * in the next "layer" of degrees of freedom. This array stores the
+ * indirect addressing of both the values and the gradient.
+ *
+ * The first table index runs through the faces of a cell, and the
+ * second runs through the pairs of the nodal degrees of freedom of the
+ * face and the derivatives, using <code>2*dofs_per_face</code> entries.
+ *
+ * @note This object is only filled in case @p element_type evaluates to
+ * @p tensor_symmetric_hermite.
+ */
+ dealii::Table<2,unsigned int> face_to_cell_index_hermite;
+
/**
* Check whether we have symmetries in the shape values. In that case,
* also fill the shape_???_eo fields.
n_q_points (0),
dofs_per_cell (0),
n_q_points_face (0),
- dofs_per_face (0)
+ dofs_per_face (0),
+ nodal_at_cell_boundaries (false)
{
reinit (quad, fe_in, base_element_number);
}
n_q_points (0),
dofs_per_cell (0),
n_q_points_face (0),
- dofs_per_face (0)
+ dofs_per_face (0),
+ nodal_at_cell_boundaries (false)
{}
fe_degree = fe->degree;
n_q_points_1d = quad.size();
- const unsigned int n_dofs_1d = fe_degree+1,
- n_q_points_1d = quad.size();
+ const unsigned int n_dofs_1d = std::min(fe->dofs_per_cell, fe_degree+1);
// renumber (this is necessary for FE_Q, for example, since there the
// vertex DoFs come first, which is incompatible with the lexicographic
scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse();
element_type = tensor_symmetric_plus_dg0;
}
+ else if (fe->dofs_per_cell == 0)
+ {
+ // FE_Nothing case -> nothing to do here
+ }
else
Assert(false, ExcNotImplemented());
// invert numbering again. Need to do it manually because we might
// have undefined blocks
- lexicographic_numbering.resize(fe_in.element_multiplicity(base_element_number)*fe->dofs_per_cell);
+ lexicographic_numbering.resize(fe_in.element_multiplicity(base_element_number)*fe->dofs_per_cell, numbers::invalid_unsigned_int);
for (unsigned int i=0; i<lexicographic.size(); ++i)
if (lexicographic[i] != numbers::invalid_unsigned_int)
{
// by reading the name, as done before r29356)
if (fe->has_support_points())
unit_point = fe->get_unit_support_points()[scalar_lexicographic[0]];
- Assert(std::fabs(fe->shape_value(scalar_lexicographic[0],
+ Assert(fe->dofs_per_cell == 0 ||
+ std::fabs(fe->shape_value(scalar_lexicographic[0],
unit_point)-1) < 1e-13,
- ExcInternalError());
+ ExcInternalError("Could not decode 1D shape functions for the "
+ "element " + fe->get_name()));
}
n_q_points = Utilities::fixed_power<dim>(n_q_points_1d);
dofs_per_cell = fe->dofs_per_cell;
n_q_points_face = dim>1?Utilities::fixed_power<dim-1>(n_q_points_1d):1;
- dofs_per_face = fe->dofs_per_face;
+ dofs_per_face = dim>1?Utilities::fixed_power<dim-1>(fe_degree+1):1;
const unsigned int array_size = n_dofs_1d*n_q_points_1d;
this->shape_gradients.resize_fast (array_size);
this->shape_values.resize_fast (array_size);
this->shape_hessians.resize_fast (array_size);
- this->face_value[0].resize(n_dofs_1d);
- this->face_gradient[0].resize(n_dofs_1d);
- this->subface_value[0].resize(array_size);
- this->face_value[1].resize(n_dofs_1d);
- this->face_gradient[1].resize(n_dofs_1d);
- this->subface_value[1].resize(array_size);
- this->shape_values_number.resize (array_size);
- this->shape_gradient_number.resize (array_size);
+ this->shape_data_on_face[0].resize(3*n_dofs_1d);
+ this->shape_data_on_face[1].resize(3*n_dofs_1d);
+ this->values_within_subface[0].resize(array_size);
+ this->values_within_subface[1].resize(array_size);
+ this->gradients_within_subface[0].resize(array_size);
+ this->gradients_within_subface[1].resize(array_size);
+ this->hessians_within_subface[0].resize(array_size);
+ this->hessians_within_subface[1].resize(array_size);
for (unsigned int i=0; i<n_dofs_1d; ++i)
{
for (unsigned int q=0; q<n_q_points_1d; ++q)
{
// fill both vectors with
- // VectorizedArray<Number>::n_array_elements
- // copies for the shape information and
- // non-vectorized fields
+ // VectorizedArray<Number>::n_array_elements copies for the
+ // shape information and non-vectorized fields
Point<dim> q_point = unit_point;
q_point[0] = quad.get_points()[q][0];
- shape_values_number[i*n_q_points_1d+q] = fe->shape_value(my_i,q_point);
- shape_gradient_number[i*n_q_points_1d+q] = fe->shape_grad (my_i,q_point)[0];
- shape_values [i*n_q_points_1d+q] =
- shape_values_number [i*n_q_points_1d+q];
- shape_gradients[i*n_q_points_1d+q] =
- shape_gradient_number[i*n_q_points_1d+q];
- shape_hessians[i*n_q_points_1d+q] =
- fe->shape_grad_grad(my_i,q_point)[0][0];
+
+ shape_values [i*n_q_points_1d+q] = fe->shape_value(my_i,q_point);
+ shape_gradients[i*n_q_points_1d+q] = fe->shape_grad(my_i,q_point)[0];
+ shape_hessians [i*n_q_points_1d+q] = fe->shape_grad_grad(my_i,q_point)[0][0];
+
+ // evaluate basis functions on the two 1D subfaces (i.e., at the
+ // positions divided by one half and shifted by one half,
+ // respectively)
q_point[0] *= 0.5;
- subface_value[0][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point);
+ values_within_subface[0][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point);
+ gradients_within_subface[0][i*n_q_points_1d+q] = fe->shape_grad(my_i,q_point)[0];
+ hessians_within_subface[0][i*n_q_points_1d+q] = fe->shape_grad_grad(my_i,q_point)[0][0];
q_point[0] += 0.5;
- subface_value[1][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point);
+ values_within_subface[1][i*n_q_points_1d+q] = fe->shape_value(my_i,q_point);
+ gradients_within_subface[1][i*n_q_points_1d+q] = fe->shape_grad(my_i,q_point)[0];
+ hessians_within_subface[1][i*n_q_points_1d+q] = fe->shape_grad_grad(my_i,q_point)[0][0];
}
- Point<dim> q_point;
- this->face_value[0][i] = fe->shape_value(my_i,q_point);
- this->face_gradient[0][i] = fe->shape_grad(my_i,q_point)[0];
+
+ // evaluate basis functions on the 1D faces, i.e., in zero and one
+ Point<dim> q_point = unit_point;
+ q_point[0] = 0;
+ this->shape_data_on_face[0][i] = fe->shape_value(my_i,q_point);
+ this->shape_data_on_face[0][i+n_dofs_1d] = fe->shape_grad(my_i,q_point)[0];
+ this->shape_data_on_face[0][i+2*n_dofs_1d] = fe->shape_grad_grad(my_i,q_point)[0][0];
q_point[0] = 1;
- this->face_value[1][i] = fe->shape_value(my_i,q_point);
- this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0];
+ this->shape_data_on_face[1][i] = fe->shape_value(my_i,q_point);
+ this->shape_data_on_face[1][i+n_dofs_1d] = fe->shape_grad(my_i,q_point)[0];
+ this->shape_data_on_face[1][i+2*n_dofs_1d] = fe->shape_grad_grad(my_i,q_point)[0][0];
}
+ // get gradient and Hessian transformation matrix for the polynomial
+ // space associated with the quadrature rule (collocation space)
+ {
+ const unsigned int stride = (n_q_points_1d+1)/2;
+ shape_gradients_collocation_eo.resize(n_q_points_1d*stride);
+ shape_hessians_collocation_eo.resize(n_q_points_1d*stride);
+ FE_DGQArbitraryNodes<1> fe(quad.get_points());
+ for (unsigned int i=0; i<n_q_points_1d/2; ++i)
+ for (unsigned int q=0; q<stride; ++q)
+ {
+ shape_gradients_collocation_eo[i*stride+q] =
+ 0.5* (fe.shape_grad(i, quad.get_points()[q])[0] +
+ fe.shape_grad(i, quad.get_points()[n_q_points_1d-1-q])[0]);
+ shape_gradients_collocation_eo[(n_q_points_1d-1-i)*stride+q] =
+ 0.5* (fe.shape_grad(i, quad.get_points()[q])[0] -
+ fe.shape_grad(i, quad.get_points()[n_q_points_1d-1-q])[0]);
+ shape_hessians_collocation_eo[i*stride+q] =
+ 0.5* (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] +
+ fe.shape_grad_grad(i, quad.get_points()[n_q_points_1d-1-q])[0][0]);
+ shape_hessians_collocation_eo[(n_q_points_1d-1-i)*stride+q] =
+ 0.5* (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] -
+ fe.shape_grad_grad(i, quad.get_points()[n_q_points_1d-1-q])[0][0]);
+ }
+ if (n_q_points_1d % 2 == 1)
+ for (unsigned int q=0; q<stride; ++q)
+ {
+ shape_gradients_collocation_eo[n_q_points_1d/2*stride+q] =
+ fe.shape_grad(n_q_points_1d/2, quad.get_points()[q])[0];
+ shape_hessians_collocation_eo[n_q_points_1d/2*stride+q] =
+ fe.shape_grad_grad(n_q_points_1d/2, quad.get_points()[q])[0][0];
+ }
+ }
+
if (element_type == tensor_general &&
check_1d_shapes_symmetric(n_q_points_1d))
{
if (check_1d_shapes_collocation())
element_type = tensor_symmetric_collocation;
else
+ element_type = tensor_symmetric;
+ if (n_dofs_1d > 3 && element_type == tensor_symmetric)
{
- element_type = tensor_symmetric;
- // get gradient and Hessian transformation matrix for the
- // polynomial space associated with the quadrature rule
- // (collocation space)
- if (fe_degree+1 == n_q_points_1d)
- {
- const unsigned int stride = fe_degree/2+1;
- shape_gradients_collocation_eo.resize((fe_degree+1)*stride);
- shape_hessians_collocation_eo.resize((fe_degree+1)*stride);
- FE_DGQArbitraryNodes<1> fe(quad.get_points());
- for (unsigned int i=0; i<(fe_degree+1)/2; ++i)
- for (unsigned int q=0; q<stride; ++q)
- {
- shape_gradients_collocation_eo[i*stride+q] =
- 0.5* (fe.shape_grad(i, quad.get_points()[q])[0] +
- fe.shape_grad(i, quad.get_points()[n_q_points_1d-1-q])[0]);
- shape_gradients_collocation_eo[(fe_degree-i)*stride+q] =
- 0.5* (fe.shape_grad(i, quad.get_points()[q])[0] -
- fe.shape_grad(i, quad.get_points()[n_q_points_1d-1-q])[0]);
- shape_hessians_collocation_eo[i*stride+q] =
- 0.5* (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] +
- fe.shape_grad_grad(i, quad.get_points()[n_q_points_1d-1-q])[0][0]);
- shape_hessians_collocation_eo[(fe_degree-i)*stride+q] =
- 0.5* (fe.shape_grad_grad(i, quad.get_points()[q])[0][0] -
- fe.shape_grad_grad(i, quad.get_points()[n_q_points_1d-1-q])[0][0]);
- }
- if (fe_degree % 2 == 0)
- for (unsigned int q=0; q<stride; ++q)
- {
- shape_gradients_collocation_eo[fe_degree/2*stride+q] =
- fe.shape_grad(fe_degree/2, quad.get_points()[q])[0];
- shape_hessians_collocation_eo[fe_degree/2*stride+q] =
- fe.shape_grad_grad(fe_degree/2, quad.get_points()[q])[0][0];
- }
- }
+ // check if we are a Hermite type
+ element_type = tensor_symmetric_hermite;
+ for (unsigned int i=1; i<n_dofs_1d; ++i)
+ if (std::abs(this->shape_data_on_face[0][i][0]) > 1e-12)
+ element_type = tensor_symmetric;
+ for (unsigned int i=2; i<n_dofs_1d; ++i)
+ if (std::abs(this->shape_data_on_face[0][n_dofs_1d+i][0]) > 1e-12)
+ element_type = tensor_symmetric;
}
}
else if (element_type == tensor_symmetric_plus_dg0)
check_1d_shapes_symmetric(n_q_points_1d);
- // face information
- unsigned int n_faces = GeometryInfo<dim>::faces_per_cell;
- this->face_indices.reinit(n_faces, this->dofs_per_face);
- switch (dim)
- {
- case 3:
- {
- for (unsigned int i=0; i<this->dofs_per_face; i++)
- {
- const unsigned int jump_term =
- this->dofs_per_face*((i*n_dofs_1d)/this->dofs_per_face);
- this->face_indices(0,i) = i*n_dofs_1d;
- this->face_indices(1,i) = i*n_dofs_1d + n_dofs_1d-1;
- this->face_indices(2,i) = i%n_dofs_1d + jump_term;
- this->face_indices(3,i) = (i%n_dofs_1d + jump_term +
- (n_dofs_1d-1)*n_dofs_1d);
- this->face_indices(4,i) = i;
- this->face_indices(5,i) = (n_dofs_1d-1)*this->dofs_per_face+i;
- }
- break;
- }
- case 2:
+ nodal_at_cell_boundaries = true;
+ for (unsigned int i=1; i<n_dofs_1d; ++i)
+ if (std::abs(this->shape_data_on_face[0][i][0]) > 1e-13 ||
+ std::abs(this->shape_data_on_face[1][i-1][0]) > 1e-13)
+ nodal_at_cell_boundaries = false;
+
+ if (nodal_at_cell_boundaries == true)
{
- for (unsigned int i=0; i<this->dofs_per_face; i++)
+ face_to_cell_index_nodal.reinit(GeometryInfo<dim>::faces_per_cell,
+ dofs_per_face);
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
- this->face_indices(0,i) = n_dofs_1d*i;
- this->face_indices(1,i) = n_dofs_1d*i + n_dofs_1d-1;
- this->face_indices(2,i) = i;
- this->face_indices(3,i) = (n_dofs_1d-1)*n_dofs_1d+i;
+ const unsigned int direction = f/2;
+ const unsigned int stride = direction < dim-1 ? (fe_degree+1) : 1;
+ int shift = 1;
+ for (unsigned int d=0; d<direction; ++d)
+ shift *= fe_degree+1;
+ const unsigned int offset = (f%2)*fe_degree*shift;
+
+ if (direction == 0 || direction == dim-1)
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ face_to_cell_index_nodal(f,i) = offset + i*stride;
+ else
+ // local coordinate system on faces 2 and 3 is zx in
+ // deal.II, not xz as expected for tensor products -> adjust
+ // that here
+ for (unsigned int j=0; j<=fe_degree; ++j)
+ for (unsigned int i=0; i<=fe_degree; ++i)
+ {
+ const unsigned int ind = offset + j*dofs_per_face + i;
+ AssertIndexRange(ind, dofs_per_cell);
+ const unsigned int l = i*(fe_degree+1)+j;
+ face_to_cell_index_nodal(f,l) = ind;
+ }
}
- break;
}
- case 1:
+
+ if (element_type == tensor_symmetric_hermite)
{
- if (this->dofs_per_face>0)
+ face_to_cell_index_hermite.reinit(GeometryInfo<dim>::faces_per_cell,
+ 2*dofs_per_face);
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
- this->face_indices(0,0) = 0;
- this->face_indices(1,0) = n_dofs_1d-1;
+ const unsigned int direction = f/2;
+ const unsigned int stride = direction < dim-1 ? (fe_degree+1) : 1;
+ int shift = 1;
+ for (unsigned int d=0; d<direction; ++d)
+ shift *= fe_degree+1;
+ const unsigned int offset = (f%2)*fe_degree*shift;
+ if (f%2 == 1)
+ shift = -shift;
+
+ if (direction == 0 || direction == dim-1)
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ {
+ face_to_cell_index_hermite(f,2*i) = offset + i*stride;
+ face_to_cell_index_hermite(f,2*i+1) = offset + i*stride + shift;
+ }
+ else
+ // local coordinate system on faces 2 and 3 is zx in
+ // deal.II, not xz as expected for tensor products -> adjust
+ // that here
+ for (unsigned int j=0; j<=fe_degree; ++j)
+ for (unsigned int i=0; i<=fe_degree; ++i)
+ {
+ const unsigned int ind = offset + j*dofs_per_face + i;
+ AssertIndexRange(ind, dofs_per_cell);
+ const unsigned int l = i*(fe_degree+1)+j;
+ face_to_cell_index_hermite(f,2*l) = ind;
+ face_to_cell_index_hermite(f,2*l+1) = ind+shift;
+ }
}
- break;
- }
- default:
- Assert (false, ExcNotImplemented());
}
}
bool
ShapeInfo<Number>::check_1d_shapes_symmetric(const unsigned int n_q_points_1d)
{
+ if (dofs_per_cell == 0)
+ return false;
+
const double zero_tol =
- types_are_equal<Number,double>::value==true?1e-10:1e-7;
+ types_are_equal<Number,double>::value==true?1e-12:1e-7;
// symmetry for values
const unsigned int n_dofs_1d = fe_degree + 1;
for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i)
for (unsigned int j=0; j<n_q_points_1d; ++j)
if (std::fabs(shape_values[i*n_q_points_1d+j][0] -
shape_values[(n_dofs_1d-i)*n_q_points_1d
- -j-1][0]) > zero_tol)
+ -j-1][0]) >
+ std::max(zero_tol, zero_tol*
+ std::abs(shape_values[i*n_q_points_1d+j][0])))
return false;
// shape values should be zero at x=0.5 for all basis functions except
return false;
const double zero_tol =
- types_are_equal<Number,double>::value==true?1e-10:1e-7;
- // check: - identity operation for shape values
- // - zero diagonal at interior points for gradients
- // - gradient equal to unity at element boundary
+ types_are_equal<Number,double>::value==true?1e-12:1e-7;
+ // check: identity operation for shape values
const unsigned int n_points_1d = fe_degree+1;
for (unsigned int i=0; i<n_points_1d; ++i)
for (unsigned int j=0; j<n_points_1d; ++j)
memory += MemoryConsumption::memory_consumption(shape_hessians_eo);
memory += MemoryConsumption::memory_consumption(shape_gradients_collocation_eo);
memory += MemoryConsumption::memory_consumption(shape_hessians_collocation_eo);
- memory += face_indices.memory_consumption();
for (unsigned int i=0; i<2; ++i)
{
- memory += MemoryConsumption::memory_consumption(face_value[i]);
- memory += MemoryConsumption::memory_consumption(face_gradient[i]);
+ memory += MemoryConsumption::memory_consumption(shape_data_on_face[i]);
+ memory += MemoryConsumption::memory_consumption(values_within_subface[i]);
+ memory += MemoryConsumption::memory_consumption(gradients_within_subface[i]);
}
- memory += MemoryConsumption::memory_consumption(shape_values_number);
- memory += MemoryConsumption::memory_consumption(shape_gradient_number);
return memory;
}
projects / dealii.git / commitdiff